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ABSTRACT 


A six degree of freedom, 20-state model of two spacecraft rendezvous is 
developed, one of which was controlled and the other considered to be passively 
tumbling. Solutions that minimize a series of performance indices are obtained for the 
problem of close approach, up to the point of contact, using a direct optimal control 
method. The solution is then verified as optimal by way of an indirect method based on 
the Minimum Principle. Next, a trajectory generation method for spacecraft reorientation 
is developed, based on a quaternion construction of Inverse Dynamics in the Virtual 
Domain. This new construction enables development of an Inverse Dynamics in the 
Virtual Domain rapid-trajectory generation method that exploits the concept of 
decoupling space and time, for the problem of a spacecraft perfonning a close approach 
maneuver to a tumbling object. Finally, the advantages of the new method are 
demonstrated through simulated scenarios that employ two distinct concepts of closed- 
loop feedback. The benefits seen by Inverse Dynamics in the Virtual Domain methods 
include the rapid computational time that allows a feasible solution to be generated, 
potentially onboard a spacecraft and in closed-loop. Although the Inverse Dynamics in 
the Virtual Domain method cannot match the true optimal solution, it has several 
advantages. Rapid computational time and the ability to reshape itself when provided 
updated state variable information reinforce the overall robustness of the method for safe 
trajectory planning. The novel trajectory generation method developed is tested using 
Monte Carlo methods to demonstrate its ability to handle realistic situations with varying 
initial conditions. 
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I. INTRODUCTION 


The rendezvous problem of two spacecraft orbiting the earth has been addressed 
in numerous publications, research thrusts and desired mission capabilities. Rendezvous 
technology has also evolved with small spacecraft development, such as the National 
Aeronautics and Space Administration’s (NASA) DART (NASA 2004), Air Force 
Research Laboratory’s (AFRL) (Experimental Satellite System) XSS-10 (Davis 2003) 
and XSS-11 (AFRL 2005), Naval Research Laboratory (NRL) SUMO (Bosse 2004), and 
Defense Advanced Research Projects Agency’s (DARPA) Orbital Express (OE) 
(Kennedy 2008). In particular, the AFRL Space Vehicle Directorate at Kirtland Air 
Force Base in New Mexico developed the XSS-11 in order to exhibit the ability for a 
small satellite to autonomously plan and rendezvous with a passive and cooperative 
Resident Space Object (RSO) in low Earth orbit (LEO). XSS-10 was a simple proximity 
mission around the upper stage that it was boosted into orbit on. XSS-11 demonstrated 
the ability to change orbits and intercept other nonrelated objects. The use of micro¬ 
satellites to inspect, service, repair, and refuel larger spacecraft is a long-term goal. The 
closest the XSS-11 approached and maneuvered around another object in space was 
approximately 500 meters. In addition, DARPA’s OE Advanced Technology 
Demonstration Program validated the technology and techniques for on-orbit refueling 
and reconfiguration of two satellites. The mission, conducted in 2007, performed several 
autonomous rendezvous and capture scenarios, including component exchange and 
propellant transfer events. The existence of these programs demonstrates that there is a 
need for a robust and effective autonomous close proximity control algorithm for 
multiple small spacecraft. Still, even with the investment in the aforementioned 
missions, spacecraft proximity operations that have any dynamic tasking attributes are 
currently executed with humans in the loop, relying heavily on an extensive operations 
staff to plan and oversee the maneuvers, and performing maneuvers while in 
communication with ground stations. This approach is extremely cumbersome with 
respect to time, manpower and overall operational footprint. 
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From a historical point of view, one of the first attempts to recapture a satellite in 
orbit was in 1992 with Intelsat6 from the Space Shuttle (Broad 1992). In this case, jet 
thrusters on both the shuttle and the satellite were repeatedly fired during the three-day 
operation to bring the two spacecraft to a complex orbital rendezvous. The ground 
controllers of Intelsat had slowed the rotation of satellite to 0.67 revolutions per minute 
so it would be easier to grasp. Still, astronauts in spacesuits tried and failed twice capture 
the satellite from orbit and move it into the payload bay of the space shuttle Endeavour. 
Success was eventually realized on the third try. 

The vast majority of previous research pertaining to operations that focused on 
rendezvous with an uncontrolled rotating object has only addressed the translational 
challenges (Jacobsen, Lee, Zhu and Dubowsky 2002; Matsumoto, Dubowsky, Jacobsen, 
and Ohkami 2003a and 2003b). This also involved segmented maneuvers and only 
matching the translational position and velocity of a point mass with that of a docking 
point on a tumbling RSO. Nolet, Kong, and Miller used a simple glideslope algorithm 
for the purpose of attempting to generate trajectories for the two-dimensional (2D) case 
with hopes of implementation on the Massachusetts Institute of Technology SPHERES 
testbed (2004). Huntington (Huntington and Rao 2008b) and Singh (Singh and Hadaegh 
2001) independently researched proximity operations but with cooperative agents. 
Henshaw mentions optimal rendezvous with a tumbling object, but there is no formal 
presentation of the results using the Minimum Principle (MP) (Henshaw 2003). 
Henshaw’s method is a variation trajectory planner capable of planning large timescale 
maneuvers (on the order of 10,000+ seconds to several years), therefore there is no 
discussion of the computational times involved to obtain a solution. 

The optimal satellite reorientation problem alone is of general interest to many in 
the field of aerospace engineering. Its interest increases as we start to consider 
autonomous rendezvous and close approach mission scenarios that couple both attitude 
and translational dynamics with the desire to match both translational and angular motion 
in preparation for docking. The available literature on optimal and efficient spacecraft 
reorientation alone is extensive (Junkins and Turner 1986; Bilimoria and Wie 1993; 
Vadali and Junkins 1984). Many civilian and military space missions need to have agile 
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attitude maneuver capability. For instance, Tactical Satellite (TacSat) 3 was intended to 
demonstrate responsive delivery of information to operational users (Davis and Straight 
2006). Due to the satellite’s Low Earth Orbit, the timeline for tasking, slewing and 
disseminating data is greatly reduced. Other challenges include the fact that the tasking 
can be modified at any moment until after a short period of time with the ground station 
who is uploading the tasking and the idea that TacSat must autonomously slew to the 
target, collect and process the data and then down link the data directly to the customer 
who is not collocated with the ground station. Current real time feedback controls are not 
optimized for minimum time (Vadali and Junkins 1984). These challenges lend 
themselves to the need for the ability to rapidly generate feasible trajectories that are 
optimized for minimum time. 

The preceding examples bring about a very interesting concept: what if only a 
feasible solution is desired, with little or no regard for maneuver time (aside from a 
maximum final time) or power used, where mission success is the only criterion 
evaluated? Is the new metric for optimization based on computational time or 
robustness? Such is the case for the reorientation of a Tactical Satellite (TacSat) (Davis 
and Straight 2006). The desire is to perform the mission within certain operating 
constraints. If the satellite happens to reorient in minimum time, nothing will be gained 
(as long as all other constraints on the dynamics are met) because the target or the 
downlink may not be in view, or the elevation mask to the target may not be cleared. For 
certain missions, due to target spacing, access times and general Concept of Operations 
(CONOPS), a hurry-up and wait mentality will not benefit the mission as much as rapidly 
generating the trajectory for the satellite to be dynamically tasked at the latest possible 
instant (not orbits ahead of time). In order for the satellite to perform the maneuver, with 
potentially target locations dynamically changing up until the point the link to the 
groundstation is broken, a feasible trajectory needs to be calculated in the minimum 
amount of time possible before the first target acquisition. In essence, there is no time to 
wait for a detailed optimal solution to be computed, but the implementation of the best 
optimized solution at that time is of utmost importance. This concept will be revisited 
throughout the different subsections of this dissertation where the rapid generation of the 
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solution and robustness of always having a feasible trajectory to follow, is considered 
more important than the perfonnance index (PI), such as minimum time, minimum 
energy, etc., to be optimized. 


The entire mission of rendezvous and docking can be broken down into several 
components and subcomponents. Figure 1, modified from Fehse (2003), illustrates the 
various segments of the mission. 


Docking 

Structural Connection 


Capture 


_t_ 

Close Range Rendezvous 

Final Approach 
Approach to capture point 
Achievement of final conditions 

Closing 

Reduction of relative distance to target 
Preparation for final approach 


_t_ 

Far Range Rendezvous 

Transfer from Phasing into close 
vicinity of target, relative navigation 

t 

Phasing 

T 

Launch 

Figure 1. Breakdown of rendezvous and docking mission elements (After Fehse 2003). 

The work in this dissertation revolves around the desire to efficiently perform 
close range rendezvous maneuvers to a tumbling spacecraft for the purpose of close 
inspection or docking. The problem of docking, contact and grappling once the chaser 
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satellite is in position is currently being explored by NRL (Creamer 2007; Tasker and 
Henshaw 2008). Among others, research by McCamish (2007) addresses the problem of 
efficiently and safely taking a satellite from far range into close proximity, or closing, to 
an RSO, albeit a nontumbling one. Therefore, this research fits in the realm of close 
approach trajectory generation, or more specifically final approach, for the purpose of 
close inspection or docking with a tumbling RSO. Launch and Phasing (maneuvering to 
match orbital elements and far approach) are assumed to have already taken place. 
Specifically, adapting the tenninology set forth by Nolet (2007), this is subset of the 
tenninal phase, starting ~10m from the RSO up until docking contact has been initiated. 
In concert with previous research in this area, it is not the intent to advance navigation 
techniques for determining the attributes of a tumbling RSO, therefore it assumes that 
exact knowledge of the RSO states are available. Complementary research efforts are 
looking at the problem of detennining the states of an RSO (Jasiobedski, Greenspan and 
Roth 2001; Tsuda and Nakasuka 2003). An exhaustive list of guidance, navigation & 
control (GNC) architecture efforts related to rendezvous can be obtained in Nolet (2007). 
The specific research addressed here currently has direct application to a specific type of 
problem scenario posed by Nolet (2007) to rendezvous with a target that is “drifting and 
tumbling” as defined below: 

Target drifting or tumbling: it is able to communicate its states, but has no 
control authority on its displacements and attitude (e.g., docking with a 
fuel-depleted target with no reaction wheels). 

Targets that are fully cooperative or able to communicate and control attitude 
information would also be candidates as considered by Romano (Romano and Hall 2006; 
Romano, Friedman and Shay 2007). This research would be coupled with advanced 
navigation methods to attempt to explore scenarios where there is no information sharing 
between satellites as listed below (Nolet 2007): 

Target not cooperating: 

• the target has no control authority, and no communication is 
occurring (e.g., docking with a dead satellite); 
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• the target is actively trying to escape the chaser, and no 
communication is occurring (e.g., military applications). 

Regardless of the specific application and overall GNC architecture, the trajectory 
generation methods will assume to have state information regarding the attitude and 
position of the tumbling RSO. 

This dissertation is organized as follows: Chapter I contains a summary of past 
and current research on this topic and an outline of several concepts employed in the 
dissertation. Chapter II focuses on the 6 DoF spacecraft rendezvous problem while 
formulating the dynamic equations for matching points of interest on two separate 
spacecraft, one of which is controlled, the other of which is tumbling. The optimal 
control problem is addressed, resulting in deriving the adjoint equations and 
transversality conditions needed to verify optimality. Solutions are obtained using a 
direct method and are verified using indirect shooting methods based on the MP 
(Pontryagin, Boltyanskii, Gamkrelidze, and Mishchenko 1964). The purpose of this 
development is to provide a baseline solution to the formulated close approach problem 
that has not yet been addressed in research. Chapter III expands on the results of Chapter 
II increasing the complexity of the problem. The same approach is taken, but the 
translational actuators (thrusters) are body mounted, the inertia matrix of the tumbling 
object is no longer taken to be identity and the end constraints require not only matching 
position and velocity of the spacecraft docking points, but also their attitude and angular 
velocity. Chapter IV introduces the concept of Inverse Dynamics in the Virtual Domain 
(IDVD) in more detail as it applies to spacecraft reorientation, a key enabler for 
operations with respect to a tumbling object that has not been addressed so far. Chapter 
V develops a method for rapid trajectory generation based on IDVD and applies to the 
spacecraft rendezvous problem. Chapter VI contains analysis and simulations of the 
trajectories generated using the methods developed. This includes comparing the 
performance to simplified, 2D cases that exist in literature and evaluating the overall 
effectiveness when subjected to the tracking challenges of current realistic systems. The 
chapter concludes with ideas for implementation in closed-loop architecture that would 
exploit the calculation and robustness advantages of IDVD methods. 
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A. OPTIMAL CONTROL PROBLEM FORMULATION FOR CLOSE 

RENDEZVOUS AND DOCKING 

Previously, the optimal control for rendezvous with a tumbling object has only 
been studied in 2D (Ma, Ma and Shashikanth 2007) with very little analysis. Studying 
the optimal solution of a given problem also entails verifying its optimality. This is done 
through rigorous analysis of the states, costates and how they relate to the MP. This 
research formulates the three-dimensional (3D) problem and utilizes a direct collocation 
method (specifically a pseudospectral method) to obtain a solution, which is then 
validated using an indirect method based on the MP. Although direct collocation 
methods enable the generation of an optimal solution (while indirect methods have a very 
limited radius of convergence), they are extremely computationally intensive and require 
a significant amount of time to converge (Yakimenko, Xu and Basset 2008, Boyarko, 
Yakimenko and Romano 2009a and 2009b). Therefore, other rapid-trajectory generation 
techniques are developed, while using the optimal solution as a baseline for comparison. 

B. INVERSE DYNAMICS IN THE VIRTUAL DOMAIN 

1. The Concept of Inverse Dynamics in the Virtual Domain 

On top of having an extremely large computational time, direct collocation 
methods may also converge to sub-optimal solution if the number of nodes is not 
sufficient (Yakimenko et al. 2008; Boyarko et al. 2009a and 2009b). In addition to that, 
the optimal solution does not have an analytical representation, which may pose problems 
when trying to implement it using a feed-forward scheme of suggested control commands 
(Yakimenko et al. 2008; Boyarko et al. 2009a and 2009b). Furthermore, the placement of 
nodes by a direct collocation method that represent the solution are not flexible can cause 
problems when trying to interpolate complex control solutions that are not necessarily 
continuous (Hurni 2009). 

This research pursues another approach exploiting the general idea of the direct 
optimization methods of calculus of variations together with an inverse dynamics 
approach. In particular, polynomials are used as basis functions to generate trajectories 
that can be traversed according to a computed analytic speed profile. An abstract 


7 



argument is used that allows the trajectory to be decoupled in space and time, while 
always resulting in a trajectory that will reach the desired endpoint conditions. The 
resulting quasi-optimal trajectory solution is capable of being generated (and updated) 
rapidly because of the reduction in the number of varied parameters due to the restriction 
on the trajectory structure by specifying a polynomial basis. Although this method lacks 
some flexibility due to the predefined structure, it provides a feasible solution that 
satisfied the endpoint constraints on the trajectory at every iteration, even when the initial 
conditions change due to disturbances or delays. Specific applications include scenarios, 
where derivative conditions on beginning and ending states need to be met, such as 
tracking missions, docking missions and other missions that are not necessarily rest-to- 
rest maneuvers. It also has application in areas, where computational time and the ability 
to immediately employ the best available solution at any given time, is heavily weighted 
as part of the overall PI of the solution. 

2. Attitude Trajectory Generation 

The goal of this research is to provide a method to determine a feasible attitude 
trajectory solution that meets endpoint requirements and dynamic constraints while 
performing a good overall maneuver relatively to a given cost function. Also, the method 
should work for any terminal conditions including nonrest to nonrest maneuvers. The 
major requirement is that the method must focus on providing a rapid, potentially real 
time solution as opposed to off-line computations even if it requires some sacrifices in 
optimality. This methods developed here are coupled with methods developed in a later 
section for the translational trajectory generation, providing the complete solution needed 
for close approach and docking with a tumbling RSO. 

3. Translational Trajectory Generation 

Trajectory generation by IDVD has been used in the past for aeronautical 
applications where the only two frames that needed consideration were the inertial and 
body frames (Yakimenko et al. 2008; Yakimenko and Siegers 2009). This research 
utilizes the concept of IDVD in the orbital frame fixed at the RSO (an intermediate frame 
between the inertial and body frame) and exploiting the Clossey-Wilshire Equations of 
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relative motion. This allows constructing the problem to analytically set the position and 
rate values of the states in the specific frame of interest, ensuring a maneuver that will 
perform the desired close approach maneuver with respect to the RSO, allowing for rapid 
trajectory reshaping or reoptimizing and reducing overall relative error. The use of a 
speed factor of a specific form (having an analytical expression for the integral of the 
inverse function), also allows flexibility in setting the node points at specific time 
intervals. This new method of rapid-trajectory generation allows the computation of an 
optimized, feasible, safe trajectory for spacecraft rendezvous requiring significantly less 
computational time than current methods. 

C. SCOPE OF THE DISSERTATION 

This dissertation advances the body of knowledge with respect to guidance of a 
spacecraft to rendezvous and dock with a tumbling object. In particular, the main new 
contributions of this dissertation are: 

1. The analytical fonnulation of the optimal control problem is presented with 
respect to a spacecraft perfonning a controlled approach in 3D to a tumbling 
object for the purpose of rendezvous and docking. A baseline optimal 
solution is presented using pseudospectral methods and verified the necessary 
conditions for optimality based on the Minimum Principle coupled with an 
indirect shooting method. 

2. An analytic quaternion trajectory is formulated using a 5th and 7th order 
Bezier Polynomial. Derivatives are formulated and applied in the virtual 
domain providing an attitude trajectory generation technique that can 
automatically match specified derivative values at the endpoints while varying 
the speed across the spatial trajectory. This fonnulation is then applied for 
rapid trajectory generation, significantly decreasing the computational time. 

3. The analytic quaternion trajectory approach based on Bezier polynomials is 
coupled with a translational polynomial scheme to create a novel trajectory 
generation scheme for the close approach with a tumbling object. The method 
is extremely flexible allowing for real time reshaping of the trajectory based 
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on dynamically changing end constraints. Also, by incorporating a speed 
factor of a particular form, this method can provide continuous mapping from 
the virtual and time domains. Analysis shows that a trajectory can be 
calculated that is within -10% of the optimal trajectory, but only requiring a 
fraction (-0.5%) of the computational time. 

4. The IDVD method is incorporated in different GNC architecture schemes to 
demonstrate its flexibility. An IDVD reshaping method is created and 
implemented in a real time rapid-trajectory generator that incorporates the 
most current knowledge of the RSO states. A hybrid method, using a 
combination of recalculating and reoptimizing the solution as well as 
reshaping based on current RSO state knowledge, is developed and 
demonstrated. 
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II. FORMULATION AND ANALYSIS OF MATCHING POINTS OF 
INTEREST IN TWO-SPACECRAFT FOR OPTIMAL 

RENDEZVOUS 

From the theoretical standpoint, this chapter elaborates on the previous work by 
Ma (Ma et al. 2007), who has studied the minimum-control effort problem for a planar 
rendezvous to a tumbling object (with only three states, coordinates x and y, and heading 
angle 6), neglecting any path constraints and relative motion dynamics pertinent to 
proximity space operations. The Sakawa-Shindo algorithm was used for calculating the 
optimal control. This chapter expands the scope by taking into account the proximity 
motion dynamics, considering the full six DoF model, and detennining both the 
minimum time and the quadratic fonnulation of minimum-control (energy) effort solution 
of the rendezvous of two satellites. It also features a comparison of the solutions obtained 
using one of the prominent direct collocation methods with the truly optimal solutions 
obtained using the MP. Another section of this dissertation takes an even further step and 
considers three different Pis, adding additional constraints to match tenninal attitude and 
angular dynamics, along with position and velocity. 

A. TWO-SPACECRAFT RENDEZVOUS MODELING AND OPTIMIZATION 

PROBLEM FORMULATION 

This section develops a model of target-chaser rendezvous. Figures 2 and 3 show 
a graphical representation of the problem. The center of the orbit frame is fixed to the 
center of mass (CM) of the tumbling RSO. The x axis points toward the Zenith. The y 
axis lies along the velocity vector of the RSO (assuming circular orbit), and the z axis lies 
along the orbit normal of the RSO. We start from the arbitrary relative position (Figure 2) 
and attempt to bring two spacecraft together for docking (Figure 3). 
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Figure 2. Depiction of the spacecraft rendezvous problem. 



Figure 3. The desired final state of two spacecraft system. 
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The set of notations to be used in this model (some of which are shown on Figures 
2 and 3) include: 


m 

df, i = x,y,z 


df, i = x,y,z 


x,y,z 


x,y,z 


o)f, i = x,y,z 


or , i = x,y,z 


o c 

co t , i = x 9 y 9 z 


O T 

co t , I =x 9 y 9 z 


qf, i = 1,2,3,4 


the mass of the chaser spacecraft; 

the position vector of the target docking point with respect to the 
target CM expressed in the target body frame; 

the position vector of the chaser docking point with respect to the 
chaser CM expressed in the chaser body frame; 

the x, y, and z position of a chaser spacecraft CM in the 
Hill’s coordinate frame; 

the x, y, and z component of the velocity of chaser spacecraft 
center of mass in the Hill’s coordinate frame; 

the x, y, and z component of the angular velocity of the chaser 

spacecraft with respect to the inertial frame, expressed in the 
chaser spacecraft principal body coordinate frame; 

the x, y, and z component of the angular velocity of the target 
spacecraft with respect to the inertial frame, expressed in the target 
spacecraft principal body coordinate frame; 

the x, y, and z component of the angular velocity of the chaser 

spacecraft with respect to the orbital frame, expressed in the chaser 
spacecraft principal body coordinate frame; 

the x, y, and z component of the angular velocity of the target 

spacecraft with respect to the orbital frame, expressed in the target 
spacecraft principal body coordinate frame; 

the components of the quaternion, representing rotation from the 
orbital to chaser body frame; 
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q ], i = 1,2,3,4 the components of the quaternion, representing rotation from the 

orbital to target body frame; 

Q the angular rate of the orbital frame with respect to the inertial 

frame. 

Using these notations, the dynamics of the two systems can now be described as 
follows (Figure 2). The translational kinematics and dynamics of a chaser spacecraft in 
the orbit frame centered at the target vehicle are presented by Hill’s equations (Vallado 
and McClain 2007): 


x = —(2Qv + 3 Q 2 x + f x ) 
m 

y = -(-2Qx + / v ) (1) 

m 

z = 1(-Q 2 z + /.) 
m 


where f x , f and f z are the applied forces (controls) expressed in Hill’s frame. 

The rotational dynamics of the chaser, described by the vector equation 
(Greenwood 1987; Wie 1988) 

I c d) c +(o c xI c co c =T (2) 


expands into the scalar quantities: 


ox 
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(7 3 C 3-/ 1 » Z C +7; C 
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In Equations (2), (3) I = diag([I^ , /£, /5 3 ]) is the inertia matrix along the 
principal axes, (o c = [m c x ,a> c y ,(D C _ ] T , and T = [T X ,T ,T_ ] T is a vector torques (symbol T 
denotes transposition). 

Similarly, for the target (we have no control of) the rotational dynamics is given 
by the vector equation 

lV+w r xlV=0 (4) 

or its scalar form 



Defining °R as the rotation matrix to convert from the body frame of the target 

{T} to the orbital frame {0} and ° c R to convert from chaser body frame {C} to the 
orbital frame {O}, we can define the angular velocity of each object (a = {T,C }) with 
respect to the orbital frame expressed in the orbital frame: 

<1 ro _ 

R (D* - 0 (6) 

a* |_Q 



Rotation matrices „R are constructed using components of the corresponding 
quaternion as follows (Greenwood 1987; Wie 1988): 

+q‘ l -q?-q? 2 (?,“??-«,“«,“) 2 

°R. 2 (q-q° + q-q;) qf-qf+qf-qf 2 (q°q° ~q’q°) V) 

2 2 (qiql+q°q’) qf-qf-qf+qf_ 
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Quaternions are a parameterization for spacecraft attitude transfonnations based 
on the fact that the general displacement of a rigid body with one fixed point is a rotation 
about a fixed axis: the eigenaxis. The definition of the quaternion is: 

/?, sin(cr / 2) 

A sinter / 2) 

q= (8) 

p 3 sin(cr/2) 
cos(cr/2) 


where cr is the scalar value of the rotational displacement about the eigenaxis and p is 
unit vector that describes the direction of the eigenaxis. The parameters of quaternions 
are propagated according to the known relation (Greenwood 1987; Wie 1988): 
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( 9 ) 


Equations (1), (3), (5) and (9) define a 20-state system of differential equations governing 
rendezvous dynamics. Combined into the state vector x these states are: 


x y, z, x, y, z, co x , to , co, , co x , co , co, ,q^ , A ? A ■ 7i ■ ■ A ■ 7. ? 7i I 


c c c c 


( 10 ) 


The governing dynamics assume six normalized controls that can be used by the 
chaser to achieve the rendezvous conditions: 


u 


/. f, f. 


f f f 

J xmax J vmax J z i 


T T 

y max z r 


( 11 ) 


For simplicity, we will further assume /) max = \m /s 2 , T. mm = \Nm , i = x,y,z. 
These controls are three normalized components of a translational force acting on the 
chaser f t , i = x,y,z, expressed in the Hill’s coordinate frame, and three normalized 
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components of a torque allowing to change chaser’s attitude, T., i = x,y,z , expressed in 
the chaser’s body frame.) All six controls are bounded: -1 < u < 1. 


Using these controls, we would like to bring the two spacecraft from some initial 
conditions, given by 20 initial values of states x j (t 0 ), i = 1,..., 20, to docking-enabling 
conditions described by matching chaser’s and target’s docking station final positions and 
velocity vectors as described by the following equations. These conditions for positions 
and velocities, expressed in the matrix form are as follows: 
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( 12 ) 


While transitioning to docking-enabling conditions, we would like to minimize and 
compare two different performance indices: 

‘f 

J = \f 0 dt (13) 

*0 


with 


/o=l 


(14) 


for minimum time and 


f 0 — + u 7 ^ + U-^ + -I(15) 

for minimum quadratic-control (or more commonly referred to as minimum energy) 
expenditure (Henshaw 2003). An equivalent formulation to the minimum time cost can 

also be expressed as J = tf, which employs an endpoint cost (Mayer) as opposed to a 
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running cost (Lagrange). A brief statement about the cost function (or PI) is in order. 
The specific cost function for minimum energy was chosen to include all actuators. The 
fact is, it can contain any subset of the controls or states. The choice for this research was 
to keep the cost function as generic as possible (including all actuator in the cost), but one 
can easily think of situations, where the cost function might be truncated to only include 
translational actuators (the case of using Control moment Gyroscopes (CMG) or Reaction 
Wheels (RW) instead of thrusters for attitude control). The specific cost function 
implemented would depend on the details of the mission, system and the discretion of the 
user, regardless of the method applied to obtain a solution. 

B. SYNTHESIS OF THE OPTIMAL CONTROL USING MINIMUM 

PRINCIPLE 

1. Formulation of the Optimal Control Problem 

We start from the general formulation for the Hamiltonian of the system with the 
state vector x, Equation (10), controls vector u, Equation (11), and running cost of 
Equations (13)-(15): 

H (k, x, u) := f 0 + (k, x) (16) 

where operator (...) on the right-hand side denotes a scalar product of two vectors, and 
k e 91 Wx is a costate vector which differential equations are to be defined later in this 
section. 

For the specific system of Equations (1), (3), (5) and (9), with the running cost 
from Equation (14) (minimum time problem) the Hamiltonian can be written with respect 
to the state vector x defined in Equation (10) as: 
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H(k,x, u) := 1 + ^X4 + A 2 x 5 + A^x 6 +— (A 4 (2Qjc 5 + 3 Q 4 +Mj) + /1 5 (-2Q.x 4 +u 2 ) + A 6 (-Q 2 x i + m 3 )) 

m v ' 


+ 


(l c -I c ) 

\ ^ 22 J 33 / 


x 8 x 9 + u 4 


A-J + 


(40 


x 7 x 9 + u 5 


A + 


\/ll 1 22) 


x 8 x 7 + u 6 


Ay 


22 


33 


(l T -I T )x X (l T -I T )x X ( f-f ) 

1/22 J 33 O n 12 1 V 33 -'ll O 10 12 1 V 11 22 / 


Ao + 


Al +' 


IXnXio 


A2 


22 


'33 


+4 |((*9 “(4 -4 -4 4 +4) Q )- X I4 "(*8 -2(Vl5 -*13*16.) Q )*15 + ( X 1 ~ 2 ( X U X 15 + *14*16 ) Q ) *16 ) 

+4t|(-(*9 "(4 -4 -4 + 4) Q )*13 +( X 8 -2(Vl5 -*13*16.) Q )*16 +4? “ 2 (x i3 * 15 + * 14 * 16 ) O) * 15 ) 
+As ^(“( x v -2(x 13 x 15 +x 14 x 16 )n)x 14 +(x 8 -2(x 14 x 15 -x 13 x 16 )Q)x 13 +(x 9 -(4 “4 -*i4 +4 ) q )*i 6 ) 

+4 ^(~(*7 -2(x 13 x 15 +x 14 x 16 )Q)x 13 -(x 8 -2(x 14 x 15 -x 13 x 16 )Q)x 14 -(x 9 -(x, 2 6 -xf 3 ~4 +4 ) q )*i 5 ) 
+A7 — (*20 ~ *17 *18 + *19 jx l8 — 4n — 2(x 18 x 19 — x 17 x 20 ) Q) x 19 + (*10 — 2 (.y 17 x 19 +x 18 x 20 )Q^x 20 j 

"'"As 2 ^( — ("*12 — (*20 *17 *is x i9 )^2^ *17 ~*"(*n ~2(A 18 x 19 — x 17 x 20 ) Q) x 20 + (*io — 2^Xj 7 X 19 +x 18 x 20 )Q)x 19 

"*"A9 (-^lO — 2 (x i7 X 19 +X 18 X 20 )Q)x i8 + (*11 — 2(x 18 Xj 9 — X 17 X 20 )Q)x i7 + (A [2 — (-^20 *17 — *18 "'" *19 ) ^ 2 ^ *20 

"*"Ao (*io — 2(x 17 x 19 +x 18 x 20 )Q)x 17 -(x n —2 (x 18 x 19 — x 17 x 20 )Q)x 18 — ^x 12 — (A 3 o —x xi — x 18 +x 19 jojx 19 j 
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(the possibility of a singular control, when A i+3 (t) = 0 , is considered in Chapter II.B.3). 

Likewise, developing the Hamiltonian for the minimum quadratic-control cost 
function, based on Equation (15) the part of the Hamiltonian that depends on the controls 
is: 


H (k, x, u) — (u^ + U -2 + Wj 2 + u 4 + w 5 + w 6 j 
+—(A 4 u, +A 5 u 2 +A 6 u 3 ) + -^A 7 u 4 +-^A s u 5 +-^A 9 u 6 


Since the controls enter Hamiltonian nonlinearly, the optimal control is not the 
bang-bang control anymore. To be more specific, the resulting optimal control that 
minimizes the Hamiltonian for the minimum quadratic-control cost function Equation 
(15) (minimum energy) is as follows: 
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u , = < 
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- 1 , 


i/: = < 
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'i+3 


I kk 

- 1 , 


X ±!<_1 

III 

_ 1 <A± 1<1 for/= 1 , 2 ,3, 
m 


X 

i 

X 


t±i> 1 
m 


( 21 ) 


/ 


i±3-<-l 


tot 


c ’ 


x., 

-1 < — 7 X < 1 for / = 4,5, 6 , where k = i- 3 




X 


?+3 


‘tt 


>1 


Now, the differential equations for costates will be the same for both optimization 
problems and are obtained via: 


i = 


'dH 

v < 3 x 


( 22 ) 


The first six adjoint equations are given by: 

X = -3X 4 Q 2 m _1 , X = 0, Aj = A 6 Q 2 m ~ l , 

X 4 = -X + 2Aflm { , X 5 = —X - 2X 4 Qm ', X 6 = -X (23) 

The next six adjoint equations, corresponding to the states 7 through 12, take the form of: 


X 7 = IU J* 4i*9 + /22 Vg ~X 3 

122 ^ 


33 


X 4 ^ + X 5 2 ^- Y 14 + X6 2 ^"13 


(24) 


The remaining eight adjoint equations, corresponding to the states 13 through 20, are of 
the form: 
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- 2x 13 Q - 2x 16 “Q - 2x 15 ~Q 


K,^\(2x n ax u -2x ti ax ti -2x t ^ K ) 

—/i] 4 — ( —(x'y — p 16 — x ]3 —x 14 + x l5 )flj 
-4 ^-(2x 15 Qx 14 + (*s -2(^15 -■*13*16) Q ) + 2 xi 6 Qx 13 + 2 x, 3 Qx i6 ) 
-^6^(2^15^13 ~( x 7 -2(^13^15 + x 14 x 16 )Q)-2x 16 Ox 14 -2x 13 ax 15 ) 


2-14 “ A3 r\ (2^14 f2 + (x 9 ( 


_| — Y“ _V 2 + V 2 

\ a 16 a 13 a 14 -r a 15 


+ O + 2x^ £2 


1 


-3 14 -(2x„nx M -2x 16 a Il5 -2x 1! ar 1( ) 


■As _ 2(w s + Vie)' Q ) + 2x 16 Qx 14 -2x 15 Qx 13 +2x 14 Qx 16 ) 

■A 6 ^(2xi 6 ^ 13 -(x 8 -2(x 14 x 15 -x 13 x 16 )Q) + 2x 15 Qx 14 -2x 14 Ox 15 ) 


As - A3 2 ( 2x l5 Qx l4 (xg 2(x 14 x 15 x 13 x 16 ) O) + 2x 15 Qx 14 2x 13 Ox 16 j 

-^ 4 ^-(2x 15 Qx 13 -2x 16 Qx ]4 +(x 7 -2(x 13 x 15 +x 14 x 16 )Q)-2x 13 Ox 15 ) 


1 


-2 ls -(2x 13 Qx, 4 -2x 14 nx, J + 2x ls nx, 6 ) 


- Aj 6 —|+2x 13 ‘Q + 2x 14 "Q — (x 9 — ( 


1 1 2 2 2.2 

-U Q -|x 16 -x 13 -x 14 +x 15 


)q) + 2x 15 2 Q 


As = -A3 ^(- 2 xi 6 Oxi 4 -2x 13 Ox 15 +(x 7 -2(x 13 x 15 +x 14 x 16 )Q)-2x 16 Qr 14 ) 
-A 4 |(2x 16 ax 13+ 2x 16 ax 13 +(x 8 -2(x 14 x 15 -x 13 x 16 )Q)-2x 14 Qx 15 ) 
-A l5 — |2 x 13 2 Q + 2x 14 2 Q-(x 9 -(xf 6 -x 2 3 -x, 2 4 +Xj 2 5 )Qj + 2x 16 2 Qj 


3 ll2 (2x„fix- 14 2x 13 nx, 4 +2x, s nx 16 ) 


(25) 


(The adjoint equations for \ 7 /i^ can be derived by swapping the associated chaser 
states and costates with the RSO variables (ex. A 3 —» A 7 ,x 13 —> x 17 ,Tj 3 —> Tj 7 ,etc.)). 
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For this problem formulation, it should also be noted that the third and the sixth 
costate equations are decoupled. This feature will be addressed further. The next step is 
to examine the transversality conditions, which with account to the following coupled 
tenninal variations, followed directly from Equation (12). The six scalar equations for 
the components of the vector e = [e 1 ,...,e 6 ]' represent the matching of position and 
velocity of the docking points in the three components of the orbital coordinate frame at 
the time t f . If positions of the docking points are matched in the orbital frame, then 

matching velocity in the orbital frame will lead to matching velocity in the inertial frame 
as well, since the transport theorem will have the same effect for coincident points. 

Let variable E represent an endpoint cost that is set to zero for both the minimum 
time and minimum-control scenarios. Introducing an auxiliary function, E , several more 
necessary conditions for optimal control (i.e., transversality conditions) can be taken into 
account as (Yan, Fahroo and Ross 2002): 

E = £ + (26) 

(=1 


f 

v 


— \T 

SET 

dx 


I t=t 


f 



(27) 


The first six equations result in expressions that only contain the parameters u / , 

i = 1.6 : 

Ti( tj -) = —o l , A 2 (tj) = — o 2 , Ay(tf) = —Vy , 

4, iff ) = ~ U A » A Vf ) = ~°5 = ~V 6 ( 28 ) 

These six equations are then properly substituted into the remaining equations to give the 
14 equations for the end conditions, described as £ = [£),...,£) 4 ]', that must be satisfied, 
along with the six docking states for position and velocity given by Equation (12) at t f , 
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and the twenty known initial conditions. For example, for the value at tf of the 7th costate 
in relation to the first transversality condition is: 


■^7 (f/0 (20i 3 *i4 +Xi5Xi 6 )x 23 2 (x i3 X 15 X \4 X 16) X 22 ) 




+ (2(Xi 3 Xi 5 - 

+ ((4 +xf 3 


X \A X \b) X 2\ (v, 6 + X |3 X 14 X| 5 )X 2 3 )i 

*^14 -^15 )"^22 — 2 ( X 13 X 14 •^ 15 ^'! 6 )-*21 ) 


O, 


(29) 


Then substituting from Equation (28) leads to: 

Cl (f/O — 'C(tf) — (2(x 13 X 14 + _ 2(x 13 X 15 — X 14 Xj 6 )x 22 ) ( — 4(^/)) 

+ (2(x 13 x 15 -X 14 x 16 )x 21 - (xf 6 + x, 2 , -xf 4 -4)x 23 ) (-1 5 4)) . (30) 

+ ((4+4 -x, 2 , -X 2 )x 22 -2(x 13 x 14 +x 15 x 16 )x 21 )(-4(t / )) = 0 


The same is done for the remaining 13 equations. 

Therefore, everything is in place to solve the minimum time Minimum Principle 
problem converted to the nonlinear programming, two point boundary value problem 
(TPBVP) numerically. Specifically, given the dynamics of the system described by 
Equations (1), (3), (5), (9), and adjoint system Equations (23)-(25), with the optimal 
controls synthesized as Equation (19) or Equation (21), with the initial states x ; (t 0 ), 

i = l,...,20, we guess on the values of the initial costates /L.(f 0 ), / = 1,..., 20 and the 
maneuver time t , (21 variable parameters), in attempt to zero six final discrepancies in 

matching final position and inertial velocity of chaser’s and target’s docking stations 
(Equation (12)), together with fourteen relations resulting from the transversality 
conditions (Equation(26)) and assure that (Yan et al. 2002) 

H(t f ) = 0. (31) 


Therefore, we have 21 varied parameters and 21 conditions to satisfy. 


2. Accounting for a Possible Path Constraint 

For the close rendezvous problem, we should obviously take into account one 

more constraint. This constraint is that the CM of the chaser spacecraft must remain at a 
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distance larger than some minimum distance (a “keep out” sphere with a radius r) from 
the CM of the RSO, which is coincident with the origin of the orbit frame. This assures 
that the chaser vehicle will not pass through the target vehicle in order to reach the 
docking position. Mathematically, this can be defined as: 

h = (Xj 2 +%2 +x] )-r 2 > 0 (32) 

where x i , i = 1,2,3 are the first three elements of the state vector (9). Furthermore, while 
a trajectory is on a path constraint h = 0 , the tangential condition must also be satisfied 
(Pontryagin et al. 1964): 

h = — - x x x 4 + x 2 x 5 + x 3 x 6 = 0 (33) 

dt 

Consequently, the Hamiltonian should be augmented by another term: 

H(k, x, u) := f 0 + (k,x) + juh (34) 


where ju is a constant and 


h = 


d 2 h 
dt 2 


(35) 


(since the path constraint has to be differentiated with respect to time twice before the 
control variables appear in the expression). The value of // is dictated as follows: 


fi = 0 if h> 0 => off the constraint boundary 

ju> 0 if h = h = h = 0 => on the constraint boundary 

Therefore, upon first contact with the path constraint, the costate values and 
Hamiltonian will be discontinuous (Bryson and Ho 1975). 

3. Possibility of a Singular Control for a Minimum Time Problem 

Upon closer inspection, we find that in the z-direction translational control u 3 is 
decoupled from all other controls. From Equation (1) we have 
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x 3 = x 6 , x 6 = m '(-Q 2 x 3 + w 3 ) 


(37) 


or 


and Equation (23) yield 


x 3 + m Li x 3 = m u 3 


(38) 


Aj = m l Cl 2 / 1 6 , \ =-Aj 


(39) 


or 

X 6 + m~ l Q 2 /L 6 =0. (40) 

Now, for the minimum time problem, the optimal control is defined as in Equation (19), 
hence for n 3 we will have 


f 1 , \<0 

= < 

3 1-1, 4>o 

With account of Equation (40), we can state that 

M 3(0 = /(^ 3 (^o)^6( ? o))- 


(41) 


(42) 


Defined by the natural frequency of Equation (40) the control u 3 can switch from 

^ 3 max = 1 to u 3 mm = — 1 or vice versa every 7TyfmCl 1 seconds. Moreover, the solution 

to Equation (40) with the initial conditions defined by \(t 0 ) and \(t 0 ) = —X i (t Q ) can 
be found analytically: 

A 6 (t) = -A 3 (t 0 )m i>5 Q 1 sin(/if 05 Qt) + /l 6 (t 0 )cos(/«~ 05 Qt). (43) 

The roots of this equation are defined as 
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t = tan 


i 


m 0 - 5 q 1 . 


(44) 


Mh) y 


From the other hand, the extremals of Equation (43) are achieved at 


t* = - tan 1 


V ^6 (t 0 ) 


0 . 5 /^,- 

m £ 2 


J 


(45) 


As seen from Equations (43)-(45), the only possibility for the singular control would be 
when 


A i (t 0 ) = 0 and A 6 (t 0 ) = 0. (46) 

In this case, A 6 (t) = 0 and the optimal control cannot be defined from Equation (41), but 
requires more rigorous analysis (Bryson and Flo 1975). 

C. METHODOLOGY FOR OBTAINING A SOLUTION AND CHECKING 

THE OPTIMALITY 

This section presents the methodology for obtaining and verifying optimal 
solutions for two problems posted in Chapter II. A. Despite the fact that the structure of an 
optimal control was defined analytically (in the previous section), it is clear that there is 
no way one could solve this problem using direct shooting approach having some 
arbitrary values of varied parameters. No matter what optimization algorithm to solve the 
TPBVP would be used, the numerical solution will diverge. That is where direct methods 
of calculus of variations become useful. In what follows, this section introduces a specific 
rendezvous scenario with the concrete specific numerical values to play with. Next, it 
describes a procedure of obtaining quasi-optimal numerical solutions for each of two 
optimization problems using direct collocation methods, also kn own as pseudospectral 
methods. Finally, a methodology of using this solution, which is very close to the true 
optimal one, to address the problem using direct shooting method for the TPBVP 
formulated in Chapter II.B.l is introduced. 
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1. Defining a Rendezvous Scenario 

This specific scenario to play with initializes the chaser CM starting at a distance 
of 5 meters from the target CM with the target having an initial angular velocity of 0.25 
rad/s in both y and z body axes without losing generality. The body coordinate frames of 
each spacecraft and the orbit frame are assumed to be coincident with the inertial frame at 
the beginning of the simulation. The chaser docking point is located at [-0.25, 0, 0] T in 
the body frame while the target docking point is located [1, 0, 0] T in the target body 
frame. The initial values of mentioned and the remaining states for computer simulations 
discussed in the following sections are presented in Table 1. The scenario also assumes 
m = 1 kg , Q = 0.005 rad/s, I' =1 = I 3x3 (where I 3x3 is the identity matrix), r 2 =1.5m and 
the maximum bound on the final time of 10 seconds. The choice of an Q value was 
made arbitrarily. 


State 

Initial 
condition 
(m and m/s) 

State 

Initial 

condition 

(rad/s) 

State 

Initial 

condition 

(quaternion) 

State 

Initial 

condition 

(quaternion) 

X\ 

0 

Xl 

0 

Xl3 

0 

X\1 

0 

X2 

5 

Xg 

0 

X\4 

0 

Xlg 

0 

*3 

0 

Xg 

0 

Xl 5 

0 

X 19 

0 

X4 

0 

Xio 

0 

X\6 

1 

X20 

1 

X5 

0 

Xu 

0.25 





X 6 

0 

Xi2 

0.25 






Table 1. The initial values of the states 


2. Solving the Problem with the Gauss Pseudospectral Optimization 
Solver (GPOPS) 

The optimal control problems posted in Section A were first solved using the 
Gauss Pseudospectral Optimization Solver (GPOPS) (Rao et al. 2009). The GPOPS is an 
open source code that implements the Gauss Pseudospectral Method for solving optimal 
control problems. As many as 150 internal nodes were chosen for the solution. Usually, 
to speed up the numerical procedure, not more than about 60 nodes are being used 
(Yakimenko, Xu and Basset 2008). The initial conditions were chosen based on Table 1 
and the final conditions were based on matching position and velocity of the terminal 
point. 
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3. Verification of the GPOPS Solution with Minimum Principle 

Differential equations for the states (Equations (1), (3), (5), (9)) and costates 
(Equations (22)-(25)), developed in the previous sections, were programmed into the 
Mathworks’ Simulink blocks as shown in Figure 4. 



M 

. 

From3 

| [*_dot] - 

x_dot 

M > - * 

u 

F " ' m<.| . 

lambda 

m —• 

w 

Fromll 

[mu] * 

From12 

mu 


Figure 4. Simulink model for integrating the states and costates using the optimal 

switching function for controls. 

The Optimal Control Law block (on the left in Figure 4), implements the control 
based on the switching conditions Equation (19) or Equation (21) developed in Chapter 
II.B.l. The upper and middle blocks integrate the states and costates, respectively, while 
block on the bottom calculates the Hamiltonian. The initial conditions for the 20 states 
are taken from Table 1, while the 20 initial conditions for costates along with the final 
time t f are guessed on. 

The far right block calculates the terminal conditions, e and C ,. The outputs are 
sent to a cost function J3 constructed as follows: 

/? = ||e(t / )|| + ||i;(t / )|| + /f(t / ) 2 (47) 
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The MATLAB unconstrained optimization function fminunc (exploiting Quasi-Newton 
method) was used to solve for the TPBVP in the following fonn: 

X=fminunc(@cost_function,xO,options,lambda,time). 

This function attempts to find a minimum of a scalar function of several variables, 
cost_function, starting at an initial estimate xO. Here, cost_function calls the 
Simulink model (Figure 4) passes the vector of initial guesses 

X = (48) 

and computes the cost function p. The options parameter of fminunc is varied 
based on the fidelity of solution required. By default, the termination tolerance the 
function value and on vector of varied parameters is set to 10' 6 . 

For the case of path constraints placed on the state variables as in Equation (31), 
20 additional parameters are needed to define the costates values at the time, t d , where 

they become discontinuous (Bryson and Ho 1975). For this case, a reset is included 
along with the integrator that will reset the costates to these new parameters if the path 
constraint is contacted. In this case, the X value is augmented to include initial guesses of 
the costate values for the time t d (Bryson and Ho 1975) and is 

X = [d 1 (t 0 ),...,T 20 (t 0 ),/i 1 (Q),...,i 20 (Q),t / ] T (49) 

The vector X is only augmented with enough costate reset values as suspected 
path constraint contact points deduced from the GPOPS solution. For example, the 
GPOPS solution suggests only one contact with the path constraint, therefore, we only 
augment the X vector with one set of A(t d ) values. 

To implement the state path constraint of the form h(\) > 0, where h does not 
depend on u, a penalty function P, was associated with the violation of the constraint that 
took the form of: 
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( 50 ) 


[ 2x t x 4 + 2x 2 x 5 + 2x 3 x 6 , h< 0 
{ 0, h> 0 


with 


P = J pdt 


(51) 


and 


J3 = \\e(t f )\\ + \\C(t f )\\ + H(t f ) 2 +P . (52) 

instead of Equation (47). This increases the cost function if the vehicle is on the 
constraint boundary and not meeting the tangency conditions stated in Chapter II.B.2. 
Otherwise, if the vehicle is not on the constraint or meets the tangency conditions while 
on the constraint, there is no penalty associated with the cost. A penalty function of this 
type is appealing because it has a smooth transition from solutions where constraints are 
not violated to solutions where they are violated. Note, that if the vehicle is on the 
boundary and meets the tangency conditions, it will not cross the boundary. 

D. OBTAINING AND ANALYZING MINIMUM QUADRATIC-CONTROL 
(MINIMUM ENERGY) SOLUTION 

The results of the GPOPS solution and then their verification with the Minimum 
Principle solution are obtained as discussed in the previous section. 

1, Minimum Energy Solution with GPOPS 

For the minimum energy rendezvous scenario set in Chapter II.A, the GPOPS 
optimization package yielded the solution shown in Figure 5. This solution, minimizing 
the quadratic-control cost returned a value of J = 0.1133. Figure 6 shows the planar 
views of the solution. The overlaid sphere is centered a the target RSO and has a radius 
equal to that of the distance the docking point of the RSO is offset from its center of 
gravity. The final maneuver time is calculated to be 10 seconds, the upper bound on the 
final allowable time for this scenario (without this limit the optimal solution would yield 
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an infinite final time). It should be noted that although position and velocities of the 
docking points coincide, the orientation of the spacecraft do not match since this 
condition was not set as a constraint. In essence, a spacecraft can be translating at a 
different rate than the docking point because the docking point is offset from the 
spacecraft CM and rotating with respect to the orbital frame. 



Figure 5. Minimum energy solution (GPOPS): The 3D view of the optimal trajectory. 

A close-up of the final position is shown in the exploded view. 

The progression of the RSO docking point position carves out a circle in the plane 
perpendicular to the angular velocity vector of the RSO because the RSO is assumed to 
have an identity inertia matrix, an assumption that is relaxed in the next chapter. 
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Figure 6. Minimum energy solution (GPOPS): the planar views of the optimal 

trajectory. 

Figure 7 shows a plot of the resulting controls time histories (solid lines) as well 
as the associated costate time histories (without units) that were used to synthesize the 
optimal control based on Equation (15) (dashed lines). Figure 8 shows the time history of 
the difference in position and velocity of the chaser and RSO docking point in the XYZ 
orbital frame, e(t), as defined by Equation (12). The time history of the terminal 
conditions due to transversality, ^(t) , are shown in Figure 9. Obviously, all the values in 
Figures 8 and 9 do approach zero as they are supposed to do. Table 2 summarizes the 
results of optimization in terms of the values of varied parameters, initial value of the 
costates and Table 3 lists the terminal values of e(t), £(t) and H(t). Since for 
numerical solutions the final value of the Flamiltonian does not tell a full story, its 
complete time history is presented in Figure 10. 
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Figure 7. Minimum energy solution (GPOPS): control and associated costate history. 



Figure 8. Minimum energy solution (GPOPS): time history of discrepancies in the 

position and velocity of the chase and RSO docking points. 
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Figure 9. 


Minimum energy solution (GPOPS): history of the transversality conditions. 


Costate 

Initial Condition 

Costate 

Initial Condition 

Costate 

Initial Condition 

X\ 

0.0202713560903818 

As 

0.0235328656506421 

Al5 

-0.00357471284206 


0.0415555942048028 

A 9 

-0.058080907113755 

Al6 

1.060304309630e-10 

A 3 

0.00694445082 087804 

AlO 

-0.358237973748141 

Al7 

-0.0681335012622741 

a 4 

0.0816306658067988 

An 

0.556788661928892 

^18 

0.0053037488085563 

A 5 

0.244702539834469 

A]2 

0.410549229662363 

A 19 

0.202974770150912 

)i 6 

0.0166838774111 

Al3 

1.41849853483e-05 

A-20 

-3.422927080667e-10 

A? 

-0.004695462475917 

Aj4 

-2.046214161772e-05 

l f 

10 


Table 2. The initial values of costates and final time as defined by GPOPS. 


endpoint 

Resulting Value 

endpoint 

Resulting Value 

endpoint 

Resulting Value 

ei(t f ) 

2.1 e-10 

cm 

-9.9 e-09 

CM 

7.7 e-08 


-9.4 e-12 

cm 

-1.6 e-09 

Ciofi/) 

5.0 e-08 

?M 

-2.8 e-11 

cm 

1.2 e-11 

CM 

1.4 e-07 

?M 

5.1 e-11 

m 

-4.5 e-09 

Cn(t f ) 

6.5 e-08 

e 5 (f) 

1.7 e-10 

Ut f ) 

2.6 e-09 

CM 

3.3 e-07 

e 6 (G 

1.4 e-10 

cm 

6.4 e-07 

CM 

-5.8 e-07 

Cl (0 

3.8 e-07 

CM 

-6.4 e-09 

m ) 

-0.009711 


Table 3. Value of terminal and transversality conditions at the final time. 
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Figure 10. Minimum energy solution (GPOPS): history of the Hamiltonian. 

It should be noted that the achieved terminal tolerance of the order of 1 O' 7 -10' 11 
does not necessarily tell about the quality of the solution. The reason for this is that the 
parameters of the trajectory are only being computed (optimality conditions enforced) at 
150 nodes. Another important issue worth mentioning here is that it took 11,869.27 
seconds (a little over 3 hours) to produce this 150-node solution of a lOs-long trajectory 
on a 2.33GHz Dell Precision M90 desktop computer with an Intel T7600 processor and 
1Gb of RAM. The initial guess for the solution (required as an input to GPOPS) consisted 
of two terminal points, one at the initial time and one at the final time. The guess for the 
initial states corresponded to the initial conditions, while the guess for the final states 
consisted of zeros for the first 12 states, and the value of [0,0,0,1] T for the states 
corresponding to the quaternions. The guess for the control history was 0 at initial and 
final time for all controls. Less accurate solution, with the lower number of nodes, would 
obviously require less computational resources (Yakimenko, Xu and Basset 2008). For 
instance, a 25-node solution, for this case, only required 219.43 seconds (less than 4 
minutes) on the same computer. The point is that the solution seems to be feasible and 
realizable in practice (look at the smooth controls in Figure 7) but can only be used for 
off-line computations, i.e., in the open-loop guidance schemes. The jump in Hamiltonian 
value in Figure 10 and other Hamiltonian histories to follow occurs when path constraint 
from Equation (32) is enforced as in Equation (36). 

The next section addresses validation of the GPOPS solution, showing that it is 
indeed optimal. 
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2. Verification of the GPOPS Solution with Minimum Principle 

As discussed in Chapter II.C.3, the initial guesses provided by the GPOPS 
solution (Table 2) were used to run an optimization procedure exploiting the optimal 
controls synthesized using the Minimum Principle. The Quasi-Newton method based 
optimization routine employing forward shooting and integration of equations of motion 
using Bogacki-Shampine ordinary differential equation 3 (ODE3) solver with a fixed step 
of 10' . This approach results in a solution that has as many as 10,000 points (as opposed 
to just 150 nodes as in GPOPS). Of course it comes with the cost. Even with the perfect 
initial guesses for all varied parameters it takes many hours for the optimization process 
to converge (compared to several minutes with no good initial guess for GPOPS). The 
MP solution, returned a value of the PI ,7=0.11341185 (compare to J = 0.1133 for the 
GPOPS solution), cost function p =0.001768, and P = 1.646 e-06. 

As expected, solving the same problem, as in Chapter II.D.l, (a quadratic-control 
case with path constraints) produces the results, the optimal trajectory, controls and time 
histories of all parameters, which are pretty close to those produced by GPOPS. The 
control time histories and those states and costates that were not shown for the GPOCS 
solution are shown. Specifically, Figure 11 shows the controls time histories, 
discontinuous at running into a path constraint at t d =9.6425, Figures 12-14 show the 
states and corresponding costates of the chaser and target RSO. The values of the initial 
costates as suggested by MP are shown in Table 4, and the tenninal values e(t f ), L,(t f ) 

and H(t ,) are shown in Table 5. Figure 15 presents the time history of the Hamiltonian. 


37 



( 7 S/W) X i (,S/LU) ( S/UJ) 



Figure 11. Minimum energy solution (MP): Control time histories and associated costates. 
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Figure 12. 


Minimum energy solution (MP): state and costate time-histories for 
translational variables of the chaser vehicle. 
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Minimum energy solution (MP): state and costate histories for the defining 
angular parameters of the chaser vehicle. 
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Minimum energy solution (MP): state and costate time histories for the RSO. 
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Figure 14. 







































































Costate 

Initial Condition 

Costate 

Initial Condition 

Costate 

Initial Condition 

k 

0.0202713560903818 

i -8 

0.0235328656506421 

il5 

-0.00357471284206 

i -2 

0.0415555942048028 

/I 9 

-0.058080907113755 

il 6 

1.060304309630e-10 


0.00694445082 087804 

ilO 

-0.358237973748141 

il7 

-0.0681335012622741 

i-4 

0.0816306658067988 

ill 

0.556788661928892 

i-18 

0.0053037488085563 

i-5 

0.244702539834469 

il 2 

0.410549229662363 

i-19 

0.202974770150912 

^6 

0.0166838774111 

il3 

1.41849853483e-05 

i -20 

-3.422927080667e-10 

h 

-0.004695462475917 

i]4 

-2.046214161772e-05 


10 


Table 4. The initial values of costates and t/as defined by MP. 


Costate 

Resulting Value 

Costate 

Resulting Value 

Costate 

Resulting Value 

e\{tf) 

0.0001 

m 

0.0001 

cm 

0.0001 


-0.0001 

cm 

-0.0001 

Cio(h) 

-0.0001 

e 3 fi/) 

-0.0006 

cm 

-0.0006 

cm 

-0.0006 

e 4 (t f ) 

-0.0033 

cm 

-0.0033 

Cn(t f ) 

-0.0033 

esitj) 

0.0012 

cm 

0.0012 

cm 

0.0012 

e b{tf) 

0.0013 

cm 

0.0013 

cm 

0.0013 

m 

-0.0001 

cm 

-0.0001 

H(t f ) 

-0.0291 


Table 5. The value of terminal and transversality conditions at the final time. 
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Figure 15. Hamiltonian for the forward shooting minimum-control solution. 


Note, some of the numbers in Table 4 appear in the bold font. These, are the only 
numbers that are different from the GPOPS solution (Table 4). Therefore, it took hours 
and hours to correct initial values for just a few costates, primarily four of them, /In, X\ 4 , 
Ti6, and a 2 o- Obviously, it is due to the fact that because if integration the solution is quite 
sensitive to the initial values of varied parameters. Also, the accuracy of GPOPS solution 
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(Table 3) cannot possibly be matched. However, the point of the entire exercise was to 
prove that GPOPS does provide a solution that is very close to the truly optimal one, and 
it was proven. 

E. OBTAINING AND ANALYZING THE MINIMUM TIME SOLUTION 
1. Minimum Time Solution with GPOPS 

The 3D trajectory and 2D projections of the trajectories as a result of the GPOPS 
solution are shown in Figures 16 and 17, respectively. As seen, they are quite different 
from those shown in Figures 5 and 6. Also notice that, although the positions and 
velocities of the angular velocities match, the vehicle orientations do not. The resulting 
control history is shown in Figure 18. The f z control (of the translational motion in the z 

direction) turns out to be highly oscillating. The endpoint conditions, taken from the state 
constraints on the endpoint as well as the transversality conditions are shown in Figure 19 
and 20. The values of the initial costates as suggested by GPOPS are shown in Table 6 
and the terminal conditions of the boundary equations are shown in Table 7. The time 
history for the Hamiltonian is presented in Figure 21. 
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Figure 16. Minimum time solution (GPOPS): The 3D view of the optimal trajectory. A 

close-up of the final position is shown in the exploded view. 
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Figure 17. Minimum time solution (GPOPS): the 2D plots of optimal rendezvous 

trajectory. 




Figure 18. Minimum time solution (GPOPS): the control and associated costate history. 
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Figure 19. Minimum time solution (GPOPS): time history of discrepancies in the 
position and velocity of the chaser and RSO docking point. 


- 1 - a in x position 

- 2 - a in y position 

- 3 - a in z position 
4 - a in x velocity 

- 5 - A in y velocity 
• 6 - A in z velocity 
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Figure 20. Minimum time solution (GPOPS): history of transversality conditions. 
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2, 

Initial Condition 

2 , 

Initial Condition 

2 , 

Initial Condition 


0.000834092040702561 

>- 

00 

-0.00781391345181247 

2-15 

-0.0670307214682799 

^2 

0.418560914087817 

Ag 

-0.0052 9954719645793 

2-16 

-2.47014407229089 

^3 

7.46232139969827e-08 

2io 

-0.0966 0150546184 92 

2-17 

-0.365359323439127 

A 4 

-0.0 06345 8 92522 04 042 

In 

0.017598938378927 

00 

0.0400467496481145 

^5 

0.901334242652277 

2 12 

-0.719562963054045 

2-19 

-0.562241844247041 

/t.6 

1.04649495108653e-07 

2l3 

-0.042 06 82 8464 05162 

220 

-1.6978671241818 

A? 

-0.0133 948 06518 9441 

2l4 

-0.0394802339637517 

*r 

3.4237 


Table 6. The initial values of costates and as defined by GPOPS. 


endpoints 

Resulting Value 

endpoints 

Resulting Value 

endpoints 

Resulting Value 


-3.1 e-08 

cm 

3.3 e-07 

CM 

-1.9 e-06 


-3.8 e-08 

cm 

-6.9 e-08 

Cioth) 

3.8 e-06 


4.9 e-08 

CM 

-5.2 e-14 

Cn(h) 

3.2 e-07 

e 4 {t f ) 

-1.8 e-07 

cm 

1.0 e-07 

Cn(t f ) 

-5.1 e-07 

es(tf) 

5.1 e-07 

cm 

1.0 e-07 

CiM 

-2.0 e-07 

e 6 {t f ) 

-1.7 e-08 

cm 

-3.5 e-07 

CiM 

1.1 e-07 

cm 

-3.3 e-10 

CM 

9.7 e-06 

mtr) 

-0.2801 


Table 7. Value of terminal conditions at the final time as calculated by GPOPS. 



Figure 21. Minimum time solution (GPOPS): history of the Hamiltonian. 


Closer inspection of the GPOPS solution for A(, in Figure 18 (it jumps back and 
forth? around zero value) suggests a presence of a singular control. This is indirectly 
confirmed by oscillations of the Hamiltonian as well (Yakimenko, Xu and Basset 2008). 
Hence, the f z control is not optimal and infeasible. 

With 150 nodes, the computational time to arrive at the solution shown above 
(3.4237 second maneuver) was 8,929.55 seconds (~2.5 hours). The initial guess for the 
solution was the same as the one used for the minimum energy case. For comparison, 

with 25 nodes the required computational time can be brought down to 100.77 seconds 

45 

















(<2 min). However, as seen from Figure 22, the controls in this case are even less 
optimal than in Figure 18, and of course the Hamiltonian (Figure 23) looks worse than 
that of Figure 21. 



Figure 22. Minimum time solution (GPOPS): The associated costate history resulting 

from a 25 node solution. 



Figure 23. Minimum time solution (GPOPS): Hamiltonian for the associated 25-node 

solution. 

2. Verification of the GPOPS Solution with Minimum Principle 

The minimum time rendezvous problem is approached in a similar fashion of the 
minimum-control one. First, the problem is investigated using a (3 based on Equation 
(52). A major difference that arises compared to the minimum-control solution is the 
existence of a singular control in M 3 , which controls acceleration in the z orbital direction. 
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As discussed in Chapter II.B.3, a singular control exists and needs special consideration 
for numerical implementation. The optimal control structure was assumed to have the 
following form: 


w 3 


-1, 0 <t<t l 

0 , t x <t <t 2 
1 , t 2 <t <t 3 
< 0, t 2 <t <t 4 
-1, t 4 <t<t 5 
0 , t 5 <t< t 6 

1 , t 6 <t<t f 


(53) 


Accordingly, the vector of varied parameters X was augmented with the switching 
times ; 


X — W(tf • (54) 

Having instances t x ,...,t 6 as varied parameters implied that the search was made 
among multiple control profdes including traditional bang-bang control, like w 3min -u 3max 
(when t x =t 2 and /,. = t f , i = 3,...,6), or w 3max -w 3min (when t x =t 2 = 0, t 3 =t 4 and 
t 5 =t 6 = t f , i - 3,...,6) allows calculation of a control history even in the presence of a 
singular arc. 

The resulting trajectory is shown in Figures 24 and 25. The corresponding optimal 
controls profdes are presented in Figure 26. As seen, it was w 3min - 0-w 3max -w 3min profde 
that was found to be optimal for the w 3 control (i.e., t 5 =t 6 = t f ). The time histories for 
the remaining controls, u l ,u 2 ,u 4 ,u 5 ,u 6 , match those for the GPOPS solution (Figure 15). 
Figures 27-29 show the state costates time histories for the chaser and target RSO, 
respectively. The values of the initial costates as suggested by MP are shown in Table 8 
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and the error in satisfying boundary conditions are shown in Table 9 (the cost /? = 0.002, 
P=0.00124 Figure 30 presents the Flamiltonian (compare it with that of Figure 21 and 
22 ). 
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Figure 24. Minimum time solution (MP): 3D optimal rendezvous trajectory. 
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Figure 25. Minimum time solution (MP): The 2D projections of the optimal rendezvous 

trajectory. 
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Figure 26. Minimum time solution (MP): control histories and associated switching 

conditions. 
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Figure 27. 
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Minimum time solution (MP): State and costate histories for the translational 
variables of the chaser vehicle. 
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Minimum time solution (MP): state and costate histories for the defining 
angular parameters of the chaser vehicle. 
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Figure 29. 


Minimum time solution (MP): state and costate time histories for the RSO. 


A. 

Initial Condition 

T, 

Initial Condition 

T, 

Initial Condition 

Xi 

0.000813270815835508 

^8 

-0.0079043944784984 

^15 

-0.067030975848676 


0.41856021740054 

Ag 

-0.0052514744594656 

A)6 

-2.47013927846309 

h 

0 

•^10 

-0.0966015055875181 

An 

-0.365359450504476 

A4 

-0.00632550328119496 

/111 

0.017598938703864 

^18 

0.0400468901069352 

A5 

0.90134410424794 

A\2 

-0.719562962786195 

A\g 

-0.562241768137901 

^6 

0 

ll 3 

-0.042050296747438 

^20 

-1.6978671241818 

I7 

-0.0135549835134852 

/ll 4 

-0.0395088309550527 


3.432 


Table 8. The initial values of costates and tf as defined by MP. 


e, 

Resulting Value 

e, 

Resulting Value 

e, 

Resulting Value 

eifi/) 

0 .016 

m 

0 .0028 

CM 

0.0096 

eit?/) 

-0.014 

cm 

-0.0018 

Cioth) 

0.0061 

Mf) 

0.0016 

cm 

0.00013 

Cn(h) 

0.00048 


-0.017 

cm 

-0.00059 

Cn(t f ) 

-0.00079 

e 5 (t/) 

0.029 

cm 

-0.00077 

Cn(t f ) 

0.0010 

ee(t/) 

0.0080 

cm 

0.0066 

C\M 

6.6 e-05 

Ci(b) 

1.9 e-05 

CM 

-0.0063 

Hit}) 

0.15 


Table 9. Value of terminal and transversality conditions at the fin al time. 
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Figure 30. Minimum time solution (MP): Flamiltonian for the forward shooting minimum 

time solution. 

Again, the bold numbers in Table 8 indicate the differences as compared to 
GPOPS quasi-optimal solution. As opposed to the minimum-control case, when truly 
optimal solution happened to have almost the same values of the initial costates, in the 
minimum time case implying a singular control, the optimal solution involved more 
variations from the solution provided by GPOPS. Of course, as in the minimum-control 
case presented in Chapter II.D. 2, the indirect-method-based optimization was run 
recursively multiple times for several days in order to converge. Nevertheless, the errors 
in satisfying the terminal conditions (Table 9) were of several orders higher than that of 
claimed by the GPOPS solution (Table 7). The next section provides some more details 
on this issue. 

3. Propagation of the GPOPS Solution 

Once again, for the solutions provided by GPOPS, vehicles dynamics are only 
satisfied in a limited number of nodes. That is why the error in meeting all constraints is 
so negligibly small as compared to solutions provided by the methods involving 
integration of equations of motion (shooting). Flowever, if we integrate the GPOPS 
controls, we will end up with about the same (lower) accuracy as the shooting approach. 

To illustrate it, let us, for example, take the controls for the minimum time 
solution (including a “weird” control for /)) and integrate it. The result of integrating the 
equations of motion derived in Chapters II. A and control switches in Chapter II.B with a 
fixed time step of 0.0001 seconds is shown in Figures 31 and 32 (the zero-order hold of 
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the last control inputs was used). As seen, the trajectories are practically the same as in 
Figures 16 and 17, but the endpoint discrepancies, summarized in Table 10, obviously 
grew up. The endpoint conditions of the state variables are shown in Table 10. Note, that 
the endpoint conditions of the transversality conditions are not known because the 
costates are not propagated. 



Figure 31. 


Propagated trajectory using the minimum time control history from GPOPS. 
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Figure 32. Propagated trajectory using the minimum time control history from GPOPS. 


e 

Resulting Value 

e 

Resulting Value 

e 

Resulting Value 

e\{tj) 

-0.00027 

Clitf) 

N/A 

cm 

N/A 

eiitf) 

-9.3e-05 

cm 

N/A 

Cio(fr) 

N/A 

^3 (tf) 

0.00037 

cm 

N/A 

Cn{t f ) 

N/A 

e\{t]) 

-0.00092 

cm 

N/A 

Cn(tf) 

N/A 

<?5 (tf) 

-0.00035 

cm 

N/A 

cm 

N/A 

ee(tf) 

0.00026 

cm 

N/A 

cm 

N/A 

cm 

N/A 

cm 

N/A 

mtf ) 

N/A 


Table 10. Value of terminal conditions at the final time. 
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III. OPTIMIZATION OF A SPACECRAFT MANEUVER FOR 
CLOSE APPROACH AND DOCKING WITH A TUMBLING 

OBJECT 


This chapter takes the study further by considering translational actuators fixed in 
the chaser spacecraft body frame (better representing an actual system), adding the 
requirement to match attitude and angular velocity of the chaser and RSO in preparation 
for a docking procedure (along with the requirement of matching position and velocity of 
chaser and RSO docking points), and considering one more PI (minimum fuel). The 
optimal control problem posed in Chapter II is, therefore, reformulated to take into 
consideration the effects on dynamics and the extra constraints. 


A. RENDEZVOUS MODELING AND OPTIMIZATION PROBLEM 
FORMULATION 


Using these controls, we would like to bring the two spacecraft from some initial 
conditions, given by 20 initial values of states -CVo), * = 1,...,20, to a docking-enabling 
condition described by matching the chaser’s and RSO’s docking station final positions 
and velocity vectors as well as matching the orientation and angular velocity of the two 

vehicles. Defining ° T R as the rotation matrix to convert from the body frame of the target 
to the orbital frame, we can express the desired terminal conditions. These conditions of 
matching spacecraft docking positions and velocities, in the matrix form, are identical to 
the conditions shown in Equation (12). Matching orientation and angular rates leads to 
seven more boundary equations at t/. 
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If Equations (55) and (56) are satisfied, that would provide: and °( 0 <: =° to / . 

Therefore, the conditions in Equation (12) can be rewritten as: 
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(58) 


The resulting complete set of desired final conditions are [e x (t f ), e 13 (/ / )] T = 0 . 

For optimal control problem in this chapter, three different perfonnance indices 
each of the (Mayer) form 

J = \f 0 dt (59) 

h 

with the running cost 


/o=l 


(60) 


for minimum time, 
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( 61 ) 


fo ~ “( M l 2 +U 2 + U { + U A +U 5 2 +U b) 

for minimum quadratic-control expenditure (i.e., minimum energy), and 

fo = (|wi | + |w 2 1 + |w 3 | + |w 4 | + |w 5 | + |w 6 1) (62) 


for minimum-control are considered. Each case has a maximum finite value chosen for 
the final time t f but is still considered a free (varied) parameter. The next step is to 

convert this optimal control problem to a TPBVP using MP (Pontryagin et al. 1964; 
Bryson and Ho 1975) 

B. MINIMUM PRINCIPLE FORMULATION 

This section deals with the MP formulation and synthesis of the optimal control. 
Using the same state vector of the complete system as shown in Equation (10), (the 
differential equations for parameters of chaser and target quaternions are standard and 
can be found in Chapter II), and the control vector shown in Equation (11), we write 
down the Hamiltonian: 

H (k, x, u) := / 0 + (k, x) + juh (63) 

where operator (...) on the right-hand side denotes a scalar product of two vectors, and 
X e S .H v ' is a costate vector where its differential equations are to be defined later in this 
section. The value of n is a constant and is dictated by the same constraint as in Equation 
(36) when considering path constraints. 

Again, as in Chapter II.B.2, this constraint is that the CM of the chaser spacecraft 
must remain at a distance larger than some minimum distance (a “keep out” sphere with a 
radius r) from the CM of the RSO, which is coincident with the origin of the orbit frame. 
This assures that the chaser vehicle will not pass through the target vehicle in order to 
reach the docking position shown mathematically in Equation (32). 

The entire Hamiltonian is expressed, with respect to the state vector x defined by 
Equation (10), in Equation (64). 
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H (k, x, u) = 1 + \x 4 + A 2 x s + A 3 x 6 

A 4 (2Qx 5 +3Q 2 Xj +(.Xj 2 6 +Xj 2 3 -x 14 -Xj 2 5 W + 2(x I3 x 14 ~x {5 x X6 )u 2 +2(x 13 x 15 + x 14 x 16 )z^ 3 ) 


H- 

m 


+/l 5 ( 2fix 4 + 2 [x n x X4 + -X 15 ^ 16 ) u x + [x l6 x 13 +x 14 x 15 )u 2 + 2(x 14 x 15 x 13 x 16 )w 3 ) 
+A 6 (-Q 2 x 3 +2(x 13 x 15 -x 14 x 16 )wj + 2(x 14 x 15 + x 13 x 16 )w 2 + (xj 2 6 -x 2 -x\ a + x 2 )« 3 ) 


v 

I c -I c 


I c -I c I c -I c 
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l 22 


+ II2 T 33 (*n*l 2 Ho + 133 tT^ 1 (Vl2)^l + 
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-(*8 -2 (*14*15 - *13*16 ) Q ) *15 + (*7 -2 (*13*15 + *14*16 ) Q )*1, 

-(^9 -((^i 6 ) 2 -(^i 3 ) 2 -(^i 4 ) 2 + (*i; 
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Q \x 
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(64) 
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For the minimum time cost problem defined by Equation (60), the part of the 
Hamiltonian that depends on the controls, the switching function, turns out to be: 

H t ('k,x,\x) = —(XJ x +X 5 f +AJ : ) + ^Z j u 4 +-^\u 5 +-^^u 6 (65) 

m hi 22 1 33 


where 


f x = (q?+<ii 2 -qf ~qf) u x + 2 (^ 2 C -qU C 4 )« 2 + 2 (m^ +^ 4 c ) m 3 
f, = 2(q?qt +? 3 C ? 4 C )*h +(q?-q?+qi 2 -q?)^ + 2 (g 2 c g 3 c -qfq*)^ ■ (66) 

fz = 2(^1 C ^ 3 L -^2^4 ) M 1 +2(^2^3 C +^1 C ^4 ) U 2 +(^4^ “ Vl " “(FT + < 7 3 C2 ) M 3 


This results in Equation (65) written as: 


* 1 , .111 
H (k,x,u) = —(switchjWj + switch,w, +switch 3 u 3 ) + — Z l u 4 + — \u 5 + — A 9 u 6 (67) 


m 


22 


33 


where the switching functions that defined the control above are given by the relations: 

switch j = A 4 (q c 4 2 + qf 2 - qf - q \ 21 ) + A s 2 [q c x q^ + q^q^ ) + \2 (qfq% - q^q^ ) 
switch, = T 4 2 ( q?q% - q^q c 4 ) + T 5 ( q c 4 2 - qf + qf - qf ) + T 6 2 ( q^qf + q^q c 4 ) . (68) 

switch, = T 4 2 (q\q 3 + q c 2 q c 4 ) + f 2 (q 2 q 4 - qfq 4 ) + A 6 (qf - qf - qf + qf ) 


Assuming there are no singular arcs, since all six controls enter the switching 
function (Hamiltonian) linearly, the optimal control that minimizes the Hamiltonian in 
the case of minimum time is the bang-bang control defined by: 


U \ 

U 4 


! 1, switchj < 0 
[-1, switchj >0 

f 1, f<0 

|-1, f>o 


1, switch, < 0 
< iu 

[-1, switch, >0 

f 1, A, < 0 

i u f, 

l-i, A>0 


1, switch, < 0 
1-1, switch, >0 
| 1 , A ,<0 

l-i, ^>0 


(69) 


Likewise, developing the Hamiltonian for the minimum quadratic-control PI based on the 
running cost (61), the switching function becomes: 
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( 70 ) 


H (k, x, u) — — (^uf + + 1/4 + w 5 + w 6 ^ 

+ — (4./1 +^5/2 +4/0 + Tc r/ ^7 w 4 + 7Z“^ m 5 +7F^ m 6 


22 


33 


Taking into account Equation ( 66 ), the resulting optimal control that minimizes the 
Hamiltonian for the minimum quadratic-control (minimum energy) PI becomes: 


u : = < 


switch; 


m 


<-l 


switch 

m 

- 1 , 


and u, = 


'+3 


- 1 , 


-iAl/ = 1,2,3 


m 

switch; 

m 


>1 


JL 


'i+3 


/ 


<-l 


(71) 


-1<^<1 / = 4,5,6, j = i -3 


A 


+3 


>1 


For the minimum fuel case, corresponding to the running cost in Equation (62), the 
switching function can be written as: 


H*(k,x,u) = (| Wj | +1 u 2 1 + |w 3 1 + |w 4 1 +1 u 5 1 + |w 6 I j 

H ( \ f2 ) "*■ ~j~c ^7 W 4 T C T c ^9 U 6 

m ' 1 1 '22 '33 


(72) 


and, therefore, the optimal control structure becomes: 
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u. = < 


and u, = 


1, 

0, 

- 1 , 

1 , 

0, 

-1, 


switch; < -m 

-m < switch; < m i = 1,2,3 
switch; > m 

4 + 3 < ' 1 a 

i= 4,5,6, 7 = 1-3 


( 73 ) 


The differential equations for costates are obtained via differentiating the Hamiltonian: 


SH 

Sx 


= -X 


(74) 


which yields: 


A l = -3A 4 fl 2 m ', A 2 =0,A 3 =A 6 Q 2 m 1 

A 4 = -A i + 2A 5 D.m l , A s = -A 2 - 2A 4 Qm 1 . (75) 

K = ~^3 


The next six adjoint equations, corresponding to costates 7-12, take the form of 
Equations (76)—(81): 


I c -I c I c -I c 

^ _ -'ll -*33 2 v i 1 22 -'ll 


7 — A s x g + AgX 8 \y OX; 6 Aj 4 Qtj 5 + \ 5 £Tx 14 + A 6 Qx 13 
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22 


‘33 


1 
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I c -I c I c - I c 

K = 33 /C 22 ^*9 + 22 /C 11 V? + 4 2 1Qx i5 - 4 2 1 ^16 - 4s 2 ' Qx i3 + 4 2 ^ 


I c -I c I c -I c 

2 _ 1 33 '22 1 , -'ll J 33 

/L 9 — T C /l ~l X 8 




C ^8 X 7 "^13 0 ^14 + 4 7 ^4 3 4 ~ ^*16 + 4 ~ ^*15 


‘22 


I -I I -I 1 1 1 1 

\0 = Jr — *^i 1 E 2 " l Jr 4*n — 4 ~f2x2o — 4 — £4 9 + 4 J'^'is 4 


‘22 


‘33 


/ -/ I -I 1111 

4 = Jr 4*12 " l Jr 4*io ~ i ~4 ~— -4 —flx 20 — 4 —Ox 17 +T 20 —Qx lg 


‘33 


(76) 

(77) 

(78) 

(79) 

(80) 
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f -I T I T -I T 1 1 1 1 

\l — ~~pr Ao*ll " l Tt _ Al-Eo — A 7 ^"^XlS As “^^17 — A 9 “^*20 Ao ^"^Xjg (^ 1 ) 


l ll 


22 


The next four equations, corresponding to costates 13-16, are of the form of Equations 
(82)—(85): 


A 3 


I 


m 


f A 4 (2 x I3 u i +2 x l4 u 2 + 2x 15 w 3 ) + /t 5 (2x 14 w 1 -2x ]3 u 2 -2x 16 w 3 ) 
+/1 6 (2x 15 Wj + 2x 16 w 2 -2x 13 w 3 ) 


“As ^-(2x 13 Ox 14 -2x 16 Qx 15 -2x 15 Qx 16 ) 


- 2x 13 2 fi - 2x 16 2 Q - 2x 15 2 Q 


■ A 4 j (“(*9 - ( x i 6 ~4 - 4 + 4) Q ) 

■4|(2x 15 Ox 14 +(x 8 -2(x 14 x 15 -x 13 x 16 )Q) + 2x 16 Ox 13 +2x 13 Ox 16 ) 
-2 16 -(2x 15 Qx 13 -(x 7 -2(x 13 x 15 +x 14 x 16 )Q)-2x 16 Qx 14 -2x 13 Qx 15 ) 


(82) 


A 4 

m 


^A 4 (-2x 14 u 1 + 2x 13 m 2 + 2 x 16 w 3 ) + /1 5 (2x 13 Wj + 2x 14 m 2 + 2x 15 u 3 )^ 
v +/1 6 (-2x 16 Wj +2x 15 m 2 -2x 14 w 3 ) 


A3 ^ (2x 14 Q + ^x 9 ( 


_| v 2 _ x 2 _ 2 2 

' a 16 a 13 •*14 ^ A 15 


j Oj + 2xj 6 ~f2 + 2 xj 5 O 


1 


- 2 , 4 -( 2 x 13 £lx 14 - 2 x, ( nx, s - 2 % £lx M ) 


■As ^( _ ( x 7 2 (x 13 x 13 +x 14 x 16 )0) + 2x 16 Qx 14 -2x 15 Qx 13 + 2x 14 Qx 16 ) 

■A 6 ^( 2 xi 6 Oxi 3 -(x 8 -2(x 14 x 15 -x 13 x 16 )Q) + 2x 15 Ox 14 -2x 14 Qx 15 ) 


(83) 
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2,5=— 

m 


f A 4 (—2x ]5 u ] -2x l( u 2 +2x 13 m 3 ) + /1 5 (2x 16 m 1 -2x 15 w 2 +2x 14 u 3 ) 


y+A 6 (2x l3 u 1 + 2x 14 w 2 + 2x 15 m 3 ) 

— ^i3 ^2x 15 Qx 14 — (x 8 — 2(x 14 x 15 — Xj 3 x 16 )Q^ + 2xj 5 Qx 14 — 2x 13 Qx 16 j 

-2, 4 ^(2x 15 Qx 13 - 2x 16 Qx I 4 + (x 7 -2(x 13 x 15 + x 14 x 16 )Q)-2x 13 Qx 15 ) 

-A,, l (2*„n*, 4 -2* 14 n* 1J + 2* ls n*, s ) 

— — |+2x 13 "Q + 2x 14 ~Q — ^x 9 — ^Xjg — x 13 — x 14 + Xj'j j Cl j + 2 xj 3 ~0 


( 84 ) 


A6= — 

m 


f A 4 (2x 16 i/j - 2x 15 m 2 + 2x 14 w 3 ) + A s (2x 15 Wj + 2 x 16 w 2 - 2 x 13 w 3 ) 


1 v +/1 6 (-2x 14 m 1 + 2 x 13 u 2 + 2 x 16 m 3 ) 

-4, 5 i(-2x, ( nv 14 -2x 15 n rii + (*, -2( W> + x„x, 6 )n)-2x M nx 14 ) 

-2 14 |(2x 16 Oxi 3 + 2x 16 Qx 13 +(x 8 -2(x, 4 x 15 -x 13 x 16 )Q)-2x 14 Ox 15 ) 
-2j 5 — |2x 13 2 Q + 2x 14 2 Q-(x 9 -(xf 6 -Xj 2 -x, 2 4 + x 1 2 5 )q)+ 2x 16 2 q| 


1 


-4 1( -(2x 1J nx 14 -2, 1J nx 14+ 2x 1! n, 16 ) 


(85) 


and, the final four adjoint equations for costates 17-20, take the form of Equations (86)- 
(89): 


2 17 = -2 I7 ^-(2x 17 Qx 18 -2x 20 Qx 19 -2x 19 Qx 20 ) 


" 2,8 — (' Y 12 _ (^20 ~ X n "^18 +-^,9)^2j 

' 2,9 — (2x 19 f2x 18 +(x n — 2(x lg x 19 — x 17 x 2 
'2 20 t(2x 19 Qx 17 — (x 10 — 2(x 17 x 19 +x 18 x. 


■ 2x 17 2 Q - 2x 20 2 Q - 2x 19 2 fi j 
) Q) + 2x 20 Qx 17 + 2x 17 Qx 20 ) 
)Q) -2x 20 Qx 18 - 2x 17 Qx 19 ) 


( 86 ) 
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Ag — —2p —|2xp O + ^Xp — ^x^o — Xp — Xp + Xp ^Oj + 2 x 20 O + 2xp Oj 


'As ^(2xi 7 axi 8 -2 x 20 Oxi 9 -2xi 9 Qx 20 ) 


■ As ^(-(*10 - 2 (*17*19 + *18*20) Q ) + 2*20^*18 ~2x [9 Qx 17 + 2x 18 Qx 20 ) 
-i 20 ^(2*20^17 -(*„ -2 (x 14 x 19 -x 17 x 20 )Q) + 2x 19 Ox 18 -2x 18 Ox 19 ) 


As = -2,7 |(-2xpQx lg -(x n -2(x 18 x 19 -x 17 x 20 )Q) + 2x 19 Ox 18 -2x 17 Ox 20 ) 
-2 1 s^( 2 *i 9 ^* 17 -2x 20 Qx 18 +(x 10 -2(x 17 x 19 +x 18 x 20 )Q)-2x 17 Qx 19 ) 


- \ 9 -^-(2x 17 Qx 18 -2x 18 Qx 17 + 2x 19 Qx 20 ) 


- 4o |(+ 2 *i7 2 ^ + 2x 18 2 Q - (x 12 - ( 


2 _ 2 _ 2 2 

I x 20 x 17 X 18 -I- x 19 


)o) + 2xp 2 Q 


Ao — A 7 2 ( 2 x 20 Qx 18 2x 17 Qx 19 + (x 10 2 (x n Xp +x 18 x 20 )Q) 2 x 20 Qx 18 j 
-As ^-(2*20^*17 + 2x 20 Qx 17 + (x n - 2 (x 18 x 19 -x 17 x 20 )Q)-2x 18 Qx 19 ) 

- As ^( 2 *i7 2 ^ + 2x 18 2 Q - (x 12 - (x 20 - Xp - xf 8 + xf 9 ) Q) + 2 x 20 2 n) 


1 


-Ao-(2 xpOxi S -2xpOx 18 +2xpQx 20 ) 


( 87 ) 


( 88 ) 


(89) 


The form shown in Equations (75)-(89) is different than presented in the previous 
chapter, again because of the body-mounted thrusters. 

The next step is to address the transversality conditions, which with account to the 
following coupled tenninal variations, followed directly from Equations (55)—(58). The 
ten scalar equations, e 1 ,e 2 ,...,e 10 , represent the matching of translational and angular 

position and velocity of the docking points in the three components of the orbital 
coordinate frame at the time t , based on Equations (57) and (58). Variable E represents an 

endpoint cost that is set to zero for both the minimum time and minimum-control 
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scenarios (Yan et al. 2002) and are equivalent to Equations (26) and (27) of Chapter II. 
The first six equations result in expressions that only contain the parameters o j , 
i=l,...,6. 


f) C 5 -^(C) -^(fo) ^3? f) fip A'sif f) ^5’ ^(c) ^6 (90) 


Realizing the specific relationships, such as: 


SE 


fix 7 


= ^7 ,f =/ 7 (VV) + y ! aIld 


fi£ 


fix 


10 




(91) 


and combining terms, seven more conditions are developed that must be satisfied, along 
with six docking states for position and velocity given by Equations (57) and (58) and 
the seven angular states based on Equations (55) and (56), atr ; .. These conditions are 

shown in Equation (92): 

(//) = ^Itf + Koif = e i4 (tf) = 9 

^56 (^/) = \tf + ^(i+3)lf + + = e (i+7)(^/) = fo f° r Z 

^7 (f/) = ^7r/ + \/+4)(f + + /(;+4)( X f/?^f/) = e (i+4)(^/) = 9? for I 


= 7,8,9 and 
= 13,. ..,16 


(92) 


that only depends on x ;/ and /. ;/ . This 

results in the remaining final conditions, C ,, along with Equations (55)—(58) that need to 
be satisfied for an optimal trajectory, £ = [<£j (t f g 7 (t f )] T = 0 . 

Other parameters for the scenario are m = 1kg, Q = 0.005 rad/s, I r =diag ([3,1,2]), 
I c = l 3x3 (where l 3x3 is an identity matrix) with inertia units of kg-m 2 and tf set to a 

maximum of 10 seconds. This scenario initializes the chaser CM starting at a distance of 
5 m from the target CM with the target having an initial angular velocity of 0.25rad/s in 
both y and z body axis. The body coordinate frames of each spacecraft and the orbit 
frame are assumed to be coincident with the inertial frame at the beginning of the 
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where /(x,,,k (/ .) is the portion of 


f 8E V 

dx 



simulation. The docking point of the chaser spacecraft is located at 0.25 m in the negative 
x body frame and the docking position of the target is located at lm in the positive x body 
frame. Therefore, the orientations of the two spacecraft will be coincident at docking. The 
resulting values of the initial states for computer simulations discussed in the following 
sections are taken from Table 1 in the previous chapter. 

C. SOLVING THE PROBLEM NUMERICALLY USING A 

PSEUDOSPECTRAL METHOD 

The optimal control problem posted in Chapter III.A was again solved using 
GPOPS (Huntington and Rao 2008). This software solves for the optimal control, state 
and costate history based on a given PI and constraints. As opposed to the 150 node 
solutions of the last chapter, 200 nodes were chosen for the solution to be consistent with 
previous research on this topic (Ma et al. 2007). 

1. Minimum Time 

The 3D trajectory and 2D projections of the trajectories are shown in Figures 33 
and 34, respectively. In these figures, the solid line with the square markers shows the 
trajectory of the docking point aboard the target RSO, the dotted-dashed line with a 
circular marker represents the chaser spacecraft’s CM, and the line with upside-down 
triangles depicts the trajectory of the docking point aboard the chaser spacecraft. The 
overlaid sphere is centered at the target RSO’s CM and has a radius equal to that of the 
distance the docking point is offset from its CM. 

The resulting control history is shown in Figure 35. Figures 35-38 show the state 
and costate histories. Figures 39 and 40 show the final difference in endpoint conditions 
approaching zero and Figure 41 shows the time history of the transversality conditions, 
Equation (92), as they approach zero. The Hamiltonian is shown in Figure 42. The initial 
values of the calculated costates are shown in Table 11. Some more details about this 
solution and required computer resources will also be provided in Chapter III.C.4 when 
comparing this solution with others. 
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Figure 33. Minimum time solution (GPOPS): 3D plot of optimal rendezvous trajectory. 
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Figure 34. Minimum time solution (GPOPS): 2D plots of optimal rendezvous trajectory. 
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Figure 35. Minimum time solution (GPOPS): control history solution. 
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Figure 36. Minimum time solution (GPOPS): state and costate histories for the defining 

translational variables of the chaser vehicle. 


68 


































































Figure 37. 
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Minimum time solution (GPOPS): state and costate histories for the defining 
angular parameters of the chaser vehicle. 



Time (s) Time (s) 


Minimum time solution (GPOPS): state and costate histories for the defining 
angular parameters of the RSO. 
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Figure 38. 
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Figure 39. Minimum time solution (GPOPS): time history of the translational endpoint 

discrepancies. 




Figure 40. Minimum time solution (GPOPS): time history of the attitude endpoint 

discrepancies. 
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Figure 41. Minimum time solution (GPOPS): history of transversality conditions. 



Figure 42. Minimum time solution (GPOPS): history of the Hamiltonian. 
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Costate 

Initial Condition 

Costate 

Initial Condition 

Costate 

Initial Condition 

A\ 

-0.0832 

A% 

-0.0300 

A\S 

1.4758 

Ai 

0.2801 

Ac, 

0.0027 

A\6 

1.9302 

Ai 

-0.0172 

^10 

-0.1521 

An 

0.7304 

A 4 

0.2251 

An 

-0.0245 

^18 

-0.9484 

A 5 

-0.4795 

An 

0.5898 

/tig 

-2.1946 

A(j 

0.0426 

^13 

-0.7500 

Ac 20 

-1.5493 

Ai 

0.1521 

A\4 

1.0625 

t f 

3.4237 


Table 11. The initial values of costates and t f as defined by GPOPS for minimum time solution. 


The important thing about this solution is that there is no singular control in 
variable u z as there was in the previous scenario described in Chapter II.E.l, where the 
translational controls were expressed in the orbital frame rather than in the body frame. 
However, the solution for T y might suggest that there may be a singular arc in the 
beginning of the maneuver. It should be noted that in case of a singular control, the 
GPOPS control output may not be feasible for onboard implementation due to it highly 
oscillating nature (Boyarko, Yakimenko and Romano 2009a). 

2. Minimum Quadratic Control (Minimum Energy) 

For the minimum energy solution, Figures 43 and 44 represent the resulting 
trajectory in both three dimensions and planar view. Figure 45 shows a plot of the 
resulting control history, as well as the associated switching function that is used for 
calculation of the control based on the running cost in Equation (61). Figures 46-48 
show the states and costates. The discrepancies in endpoint conditions are shown to 
reach zero at tf in Figures 49 and 50. The transversality conditions are shown in Figure 
51 and the Hamiltonian is shown in Figure 52. Table 12 summarizes the initial values of 
the costates that were solved for by the GPOPS solution. The PI of the solution was J = 
0.2445. 
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Figure 44. Minimum energy solution (GPOPS): 2D plots of optimal trajectory. 
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Figure 45. Minimum energy solution (GPOPS): control history with respect to the 

optimal trajectory. 




Figure 46. 


Minimum energy solution (GPOPS): state and costate histories for the 
translational variables of the chaser vehicle. 
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Figure 47. 
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Minimum energy solution (GPOPS): state and costate histories for the defining 
angular parameters of the chaser vehicle. 



Time (s) Time (s) 

Minimum energy solution (GPOPS): state and costate histories for the defining 

angular parameters of the RSO. 
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Figure 49. Minimum energy solution (GPOPS): history of the translational endpoint 

discrepancies. 
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Figure 50. Minimum energy solution (GPOPS): history of the rotational endpoint 

discrepancies. 
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Figure 51. 


Minimum energy solution (GPOPS): history of the transversality conditions. 


x 10' 4 



Figure 52. Minimum energy solution (GPOPS): history of the Ffamiltonian. 

The solution GPOPS arrives at is feasible in that it does not violate any of the 
constraints presented in the problem formulation. Furthermore, the control history 
follows the logic that was derived from Minimum Principle stated in Equation (23). 
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Costate 

Initial Condition 

Costate 

Initial Condition 

Costate 

Initial Condition 

A\ 

- 0.0335 

h 

0.0225 

Al 5 

- 0.1798 

h 

0.0602 

Ag 

0.1005 

Al6 

0.0457 

h 

- 0.0356 

2 io 

0.0460 

An 

- 0.2236 

X4 

- 0.1013 

An 

- 0.1932 

2-18 

- 0.0086 

A5 

- 0.2672 

An 

- 0.1462 

2-19 

- 0.1366 

2-6 

- 0.0714 

2|3 

0.1229 

2 20 

- 0.3664 

h 

- 0.0460 

A \4 

- 0.1501 

c 

8.8943 


Table 12. The initial values of costates and t f as defined by GPOPS for minimum energy 

(quadratic-control) solution. 


3. Minimum Absolute Control (Minimum Fuel) 

To avoid using an absolute value (which is a nondifferentiable function at the 
value zero), the running cost in Equation (62) was substituted with: 

f 0 = {U\-\-U 1 U ^ + (93) 


with the corresponding modifications in the control allocation: 


u 


r rri rri rri rri rri 

lluuulll 

x max. /max /max xmax /max zmax xmax /max /max xmax /max zmax 


(94) 


so that now 0 < u < 1 (the superscript on the control, + or -, denotes the direction in the 
body axis). This also results in the following control relations (compare with Equation 
( 21 )): 


fx (#4 (7l 9l ?3 )(**1 *b) 2 (/j #2 9l 9i, )(^2 2 (/] #3 + C / 2 9 4 ) (M 3 U 9 ) 

fy= 2 ( 9 x 9 i + 4 )(«, -«7) + (^4 C2 -^1 C2 + # - 9 } 2 )(M 2 -M 8 ) + 2(<? 2 C </' - qfq 4 ) (u 3 - w 9 ) 

/z = 2 (^ q C i -q C 2 q 4 ){u x -u 1 ) + l(q;q 4 + q[q 4 )(w 2 -u 8 ) + (<?f -^ C2 -+ #)(« 3 “«#) • ( 95 ) 

T x = (m 4 -m 10 ) 

T y =(u 5 -u n ) 
r =(u 6 -u n ) 


Analyzing the structure of the switching function, 
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H (k,x, u) ~^ j u i +—(A 4 f + A 5 f 2 + + — 

III 1, 


1 


1 


c A 1 u 4 + c A%u 5 + c A 9 u 6 

-'ll -*22 -*33 


( 96 ) 


and comparing it with that of Equation (26), we can derive the optimal control strategy as 
follows: 


and 



and 



switch,. 

< —m 

switch,. 

> —m 

4+3 < 

U :Si 

11 — 1 

1 

Al 

m 

+ 

NT 

~!jj 

switch,. 

•6 ^ 111 

switch,. 

6 > m 

4-3 

U 

* —■ 

VI 

4-3 

>45 


/ = 1,2,3 
i = 4,5,6, j = i- 3 

i = 7,8,9 


/ = 10,11,12, y =1—9 


(97) 


Here, switch,., for i = 1,2,3, are the same as in Equation (68). Hence in this case we can 
expect bang-off-bang control. 

The results are presented in Figures 53-62 and Table 13. The value of the PI was 
computed as J = 2.4851. In Figure 20, the controls m. , / = 1,...,6, acting in the positive 

direction along the corresponding body axis appear as positive values, and u j , i = 7,..., 12 
are shown with the negative sign. Note, that according to the optimal control structure of 
Equation (97) the switching occur when the switches exceed a certain value, not zero as 
for example in Equation (69). For instance, for u z the motion starts with some short 
control input in the negative direction, corresponding to u~, and ends with some short 
control input in the positive direction, corresponding to u 4 (during the remaining time the 
control u. is zero). Also, the body mounted translational actuator, u x , shows rapidly 

oscillating control values near the end of the maneuver, which is a property of a singular 
arc. The resulting rapidly oscillating control in that region would make the solution 
incredibly difficult to implement. 
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Figure 53. 


Minimum fuel solution (GPOPS): 3D plot of optimal rendezvous trajectory. 
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Figure 54. Minimum fuel solution (GPOPS): 2D plots of optimal rendezvous trajectory. 
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Figure 55. Minimum fuel solution (GPOPS): control history with respect to the optimal 

trajectory. 
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Figure 56. Minimum fuel solution (GPOPS): state and costate histories for the 

translational variables of the chaser vehicle. 
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Figure 57. 
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Minimum fuel solution (GPOPS): state and costate histories for the defining 
angular parameters of the chaser vehicle. 
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Minimum fuel solution (GPOPS): state and costate histories for the RSO. 
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Figure 58. 
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Figure 59. Minimum fuel solution (GPOPS): history of the translational endpoint 

discrepancies. 
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Figure 60. Minimum fuel solution (GPOPS): history of the rotational endpoint 

discrepancies. 
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Figure 61. Minimum fuel solution (GPOPS): history of the transversality conditions. 
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Figure 62. Minimum fuel solution (GPOPS): history of the Ftamiltonian. 


84 

















Costate 

Initial Condition 

Costate 

Initial Condition 

Costate 

Initial Condition 

Ai 

- 0.2725 

A& 

- 1.0009 

A\5 

1.6243 

Ai 

0.4399 

Ag 

1.0180 

Ak, 

0.9074 

Ai 

- 0.1960 

A\o 

0.7419 

An 

1.1494 

A4 

0.1341 

An 

- 0.3327 

An 

- 3.2616 

As 

- 1.3660 

An 

- 0.1168 

Aig 

- 4.2237 

Ae 

- 0.8983 

An 

- 1.5419 

■A 20 

- 2.4355 

Ai 

- 0.7419 

Al4 

2.6216 

C 

8.1182 


Table 13. The initial values of costates and t f as defined by GPOPS for minimum fuel solution. 


4. Solution Comparisons and Propagation 
a. Solution and Comparison 

Table 14 summarizes the results for each 200-node solution presented in 
Chapter III.C. Obviously, the scenario with the minimum time PI uses the most fuel to 
complete the maneuver. The minimum energy solution has an order of magnitude lower 
cost, as it spreads the control over the entire maneuver. It also exhibits a smoothest 
controls time histories, but maneuver takes the longest to complete. The minimum fuel 
solution results in about the same duration of the maneuver as the minimum energy 
solution, but uses about 30% less fuel. As opposed to minimum energy, it exhibits a 
bang-off-bang control structure that was expected with this type of the running cost. 



Minimum time 

Minimum-Quad- 

Control 

Minimum fuel 

Time of the maneuver 

3.5086 

8.8943 

8.1182 

Energy expenditure 

10.4471 

0.2445 

1.1587 

Fuel expenditure 

20.9419 

3.6941 

2.4863 

5jc 

Computation time (sec) 

33,370 

86,040 

84,261 


Table 14. Comparison of solutions for the three optimal control problems. 


Table 14 also presents the computational time it took to converge to a 
solution from the very same initial guess. All computations were performed on a 
2.33GHz Dell Precision M90 desktop computer with an Intel T7600 processor and 1Gb 
of RAM. As compared to the results discussed in Chapter II (Boyarko et al. 2009a), the 
required central processing unit (CPU) time is about three times longer. This is due to the 
fact that 200 nodes were used in the simulation as compared to 150 nodes to that of 
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Chapter II. As explained by Yakimenko (Yakimenko, Xu and Basset 2008), the required 
CPU time can be drastically decreased by decreasing the number of nodes, but apparently 
the lesser-node solution would probably not catch short impulses of Figure 47. 

The ultimate goal is to allow trajectory generation to be carried in real 
time onboard spacecraft. At this point, this cannot be realized using either pseudospectral 
methods (Yakimenko et al. 2008). 

b. Propagating the Solution 

The goal of this section is to check the feasibility of using the optimal 
solution in order to control the system in feed-forward mode. Consider the very first 
case, when the minimum time solution was obtained. Now that we have the controls, let 
us integrate the equation of motion derived in Chapter III.A with a time step of 0.0001 
seconds and use the calculated controls as an input. The resulting trajectory is shown in 
Figures 63 and 64. The time history of the differences in chaser and RSO docking 
position, velocity, angular position and angular rate are shown in Figure 65. Table 15 lists 
the resulting state variable deviations from the desired endpoints stated in Equations (55)- 
(58) at the final time. 

The propagated trajectory appears to be very close to that of Figures 32 
and 32, i.e., propagation of the dynamics using the minimal time control output from 
GPOPS results in a feasible trajectory that is similar to the resulting trajectory generated 
by GPOPS. This means that if the solution could be obtained in real time, with a perfect 
model of the real system, it can be passed to the control system by feed-forwarding the 
controls time histories. 
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Figure 63. Propagated trajectory using the minimum time control history from GPOPS. 
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Figure 64. Propagated trajectory using the minimum time control history. 
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Figure 65. Endpoint conditions on chaser and RSO docking position, velocity, angular 

position and angular rate. 


e 

Resulting Value 

e 

Resulting Value 

e 

Resulting Value 


-0.0003 


0.0001 

cm 

N/A 

eiUf) 

-0.0021 

e9 {tj) 

-0.0001 

cm 

N/A 

<?3 (tj) 

-0.0007 

e\o{tf) 

0.0005 

cm 

N/A 

e 4 (t f ) 

-0.0003 

eu (tj) 

0.0004 

cm 

N/A 

esitf) 

0.0015 

e n(tj) 

-0.0002 

cm 

N/A 

ee(tf) 

-0.0014 

ei3 (tf) 

- 0.0000 

cm 

N/A 


-0.0001 

cm 

N/A 

H(tf) 

N/A 


Table 15. Value of terminal conditions at the final time. 
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IV. TIME-OPTIMAL REORIENTATION OF A SPACECRAFT 
USING A DIRECT OPTIMIZATION METHOD BASED ON 
QUATERNION REPRESENTATION OF THE INVERSE 

DYNAMICS 

This chapter of the manuscript focuses on using the equations of motion to define 
the controls for a given attitude trajectory. This is the first step in defining the necessary 
equations and relations for the IDVD method for a close approach that would enable 
inspection or docking. In later sections, the methodology will be applied to the 
translational portion of the problem, but first the attitude control is considered. An initial 
trajectory is provided as an initial guess. The trajectory is then perturbed in the solution 
space by varying higher order derivatives on the endpoints until an acceptable solution is 
found. It should be noted that using inverse dynamics to optimize a rotational motion of 
a satellite has been already evaluated by other authors as well (Louembet, Cazaurang, 
Zolghardi, Charbonnel, and Pittet 2007; and Mclnnes 1998). However, the Euler angles, 
suffering from well-known kinematic singularities were used, and no attempt to decouple 
the domains of space and time were made. Furthermore, this research extends the 
situations to more realistic scenarios where nonzero rates and accelerations are present as 
in cases of target tracking and docking. 

A. SPACECRAFT MODEL AND ATTITUDE TRAJECTORY OPTIMIZATION 
PROBLEM 

The first step is to define the dynamics of the system in question, a reorienting 
spacecraft. The rotational dynamics of the spacecraft can be described by Euler’s 
rotational equations of motion. Written in the body-fixed principal axis, this results in the 
vector equation (Greenwood 1987; Wie 1988): 

Id) + oj x Ito = T, (98) 

which expands into the scalar equations: 
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go, = a 


(I22 In)°- > y°- > z T 


a X 


co„ = a = 

y y 


co„=a_ = 


(I 33 -I u )qJxO0 2 +T. 

122 

(I n -I 22 )co v 6J x + T 


33 


( 99 ) 


In turn, rotational kinematics can be described using quaternions, q = [c/, ,q 2 ,q 3 ,q 4 ]'. The 

definition of the quaternion is taken from Equation (8) where a is the scalar value of the 
rotational displacement about the eigenaxis and p is unit vector that describes the 
direction of the eigenaxis. The dynamics of the unit quaternion are described as 
(Greenwood 1987; Wie 1988): 


9 , 


0 

0} z 

~ «f 

Of 


q 2 

1 

- co z 

0 

of 

°f 

q 2 

^3 

~ 2 

°f 

- ( f 

0 

co z 

q 3 

Aa\ 



~ °f 

- 

0 

3a 


( 100 ) 


In Equations (99) and (100), I = diag([I n ,I 22 ,I 33 ]) is the inertia matrix (along the 
principal axes), to = [co x , co , co. ] T is the vector of angular velocities, and T = [T x , T y , 7) ] ' 
is the vector torques (bounded controls). 

The problem in question is to find the slew trajectory (quaternion time history) for 
a satellite subject to specific constraints that minimizes the time to complete the 
maneuver, t f . This is expressed as minimizing the PI: 


J = | dt , 

o 


(101) 


while reorienting a satellite from the initial conditions to 0 , q 0 to final conditions o> f , q 
for a system (99)-(100), subject to constraints on controls 


T < T < T 


( 102 ) 
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Bilimoria and Wie (1993) formulated this problem for a rest-to-rest maneuver and 
presented the solution using an indirect method. They showed that, in general, for the 
case of symmetric body with bounds on each torque component, it does not result in an 
eigenaxis maneuver. In addition to that, the following section presents even a more 
general solution obtained off-line to be used along with that of Bilimoria and Wie (1993) 
as a reference one for the proposed on-line solution obtained using a direct method 
exploiting inverse dynamics of Equations (99)-(100). 

Table 1 shows the different test cases examined in this chapter. Test Case 1 was 
taken directly from Bilimoria and Wie (1993), while the others were chosen to illustrate 
different scenarios of interest. 


Case 

Normalized Inertia Matrix 

Case 1 

I = diag([\,\,\]) 

Case 2 

l = diag{[3,\,2]) 

Table 16. 

Description of the test cases. 


This chapter considers two basic scenarios assuming a 0 =90° and 0=180° slew 
maneuver about the z-axis (so that q 0 =[0,0,0,l] r and q, = [O,O,sin|0,cosy0] r ) with 
zero and nonzero normalized body rates at endpoints: a) to 0 = to / = 0 3xl , b) 
to 0 = —O) ; = 75 - 1 3x |, and c) to 0 = -to, = l 3xl . For the normalized states, the constraints (5) 
take the form -l 3xl < T < l 3xl . 

B. SOLVING THE PROBLEM USING THE GAUSS PSEUDOSPECTRAL 
METHOD 

First, the problem formulated in Chapter IV.A is addressed using the same 
pseudospectral direct method as in Chapters II and III. The goal is to have some reference 
solutions for comparison to the solutions obtained through other methods, specifically 
IDVD. 

The optimal (Bilimoria and Wie 1993) solution has been matched using a 
different optimization method in previous work (Fleming 2004). 
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Figures 66-68 show the results of applying GPOPS to obtain minimum time 
solutions for a 180° slew of a satellite. Specifically, Figures 66 and 67 present time 
histories of all states and controls for the solution that involves 100 nodes, which results 
in (100-2)x 10=980 variable parameters, derived from the 98 internal node points for each 
of the 7 state and 3 control variables. Figure 68 depicts the 3D representation of the 
solution in inertial space, clearly showing that it is not an eigenaxis maneuver, with an 
inclination of the z b axis during rotation in the x b y b direction. 

This solution compares with the solution presented in Bilimoria and Wie (1993) 
fairly well. The final calculated maneuver time, t f , was found to be 3.243 seconds. 

However, it took almost two hours of CPU time to obtain this solution. Another 
observation is that because of the nature of the system described by Equation (99), the 
optimal control has a bang-bang structure as shown in Figure 67. That results in the 
maximum magnitude of the angular acceleration at the boundary points (for Case 1 
angular accelerations are simply equal to the corresponding controls). It means arriving at 
the terminal conditions with the maximum angular acceleration. Also, if we have to 
update a trajectory while the satellite performs this rotation (to accommodate possible 
disturbances and unmodeled dynamics), it would cause discontinuities of angular 
acceleration (sudden jumps in controls). Increasing the order of the system to account for 
the boundary conditions on angular accelerations will obviously cause a slight 
degradation of the PI and a further increase of the required CPU time to obtain a solution. 
Hence, although in this case GPOPS does produce a valid solution, as previously stated in 
Chapters II and III, it is not practical and currently cannot be used for online 
computations. 
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Figure 66. Case 1 (GPOPS solution): time histories of the state variables. 
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Figure 67. 


Case 1 (GPOPS solution): time history of the controlling torques. 
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Figure 68. Case 1 (GPOPS solution): the 3D representation of the solution. 

As pointed out in GPOPS documentation (Rao et al. 2009), reducing the number 
of nodes may lead to a more robust (in terms of computational time) result, therefore an 
attempt was made to obtain a solution of the same problem using a lesser number of 
nodes. These GPOPS solutions are shown in Figures 69-71. 

It turns out that for a lesser number of nodes the GPOPS converges to different 
solutions. To this end, Figure 68 shows time histories of the angular velocity components 
for the 25- and 50-node solution (involving 230 and 480 varied parameters, respectively). 
Obviously, they are different from that of the 100-node solution in Figure 64. While a 25- 
node solution is simply symmetrical with respect to the 100-node solution, as can be seen 
by comparing Figure 66 and Figure 68 showing an inclination of the z h axis during 

rotation in the -x b y b direction, and represents another equally optimal solution out of 
possible four solutions (Fleming 2004). This is due to the equivalent possibilities of 
positive or negative inclination in the z h direction coupled with possibilities of -x h y h 

and x b y b directional rotation for the 180° maneuver about a principal axis. The 
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possibility of the positive or negative inclination in the z b direction also leads to two 

equivalent solutions for similar slews less than 180°. A 50-node solution, shown in 
Figure 69, appears not to be valid (optimal) at all. 

As expected, decreasing the number of nodes leads to a substantial decrease in the 
computational time, but as shown above the method could produce a nonvalid solution. 
Also, even if it produces a valid solution the time histories for control torques, Figure 70, 
may not be trackable by the inner-loop controllers. Therefore, the solutions obtained by 
GPOPS may not be used in a real time feed-forward control scheme. 




Figure 69. Case 1 (GPOPS solution): comparison of time histories of angular velocity 
components obtained for 25 nodes (top) and 50 nodes (bottom). 
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Figure 70. Case 1 (GPOPS solution): comparison of time histories of torques, obtained 

for 25 nodes (top) and 50 nodes (bottom). 



Figure 71. Case 1 (GPOPS solution): the 3D representation of the 25-node solution. 

For Case 2, the nonsymmetrical inertia with the bounds on individual control 
torques, the solution is slightly different (Figures 72-73). The overall characteristic of 
this and other solutions involving different sets of the boundary conditions will be 
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presented in Chapter IV.E, but the general tendency is the same—it requires at least a 
hundred nodes to produce a valid and feasible off-line solution. Yet, GPOPS again 
presents a good and easy-to-use tool to produce reference trajectories that can be used for 
comparison with solutions obtained using other approaches. 




Figure 72. Case 2 (GPOPS solution): time histories of the state variables. 
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Figure 73. Case 2 (GPOPS solution): time history of the controlling torques. 
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C. APPLICATION OF INVERSE DYNAMICS IN THE VIRTUAL DOMAIN 

APPROACH WITH QUATERNION ATTITUDE REPRESENTATION 

One of the two main ideas of Inverse Dynamics in the Virtual Domain (IDVD) 
method is exploiting the differential flatness property of the equations of motion 
(Yakimenko 2000; Cowling, Yakimenko, Whidborne and Cooke 2007; Yakimenko and 
Siegers 2009). In the above problem, this relates to the fact that all the state and control 
variables can be expressed as explicit functions of the output variable or time derivatives 
of the output variable, which in this case is the quaternion itself: 

®=/i(q.q)» T^Cq^q)- (103) 

Another aspect of IDVD involves handling computations in the virtual domain allowing 
space and time decoupling. By doing so, a trajectory can be computed while also 
manipulating the speed at which that trajectory is followed. The following sections 
present a novel parameterization for the output variable, components of the quaternion, 
and develop a step-by-step computational scheme. 

1. Quaternion Parameterization 

In order to parameterize the problem, the output trajectory is approximated using 
some combination of basis functions. The standard approach would be to choose some 
combination of polynomials or trigonometric functions for the output variables 
(Yakimenko 2000; Cowling et al. 2007; Yakimenko and Siegers 2009). While this may 
be straightforward when dealing with state variables in translational space, it may 
become more challenging when dealing with expressions for orientation. 

While a quaternion may be the preferred method to express attitude because of the 
lack of singularities, choosing basis functions becomes more challenging because a 
nonlinear unit norm condition needs to be preserved across the quaternion history (Milam 
2003). Simply choosing arbitrary polynomials for the 1st three components of the 
quaternion and attempting to enforce the quaternion constraint is insufficient because it 
cannot be guaranteed the constraint will be met throughout the maneuver (specifically if 
one of the iterations on a polynomial provides a norm > 1 in any one of the components), 

so simply tacking on a penalty function for violating the constraint will not work. For 
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this reason, a specific polynomial expression for the quaternion was chosen inspired by 
the work of Kim, Kim and Shin (1995). This consists of expressing the quaternion time 
history as an exponential function containing a constant parameter multiplied by a Bezier 
polynomial of a degree n: 

q( r )=q o ri ex P(®A«( r ))’ ( 104 ) 

i =1 

where 

A»=X/U r ) < 105 > 

j=i 

and 

P in {r) = [ n \\-ry-'r 1 . (106) 

A complete list of quaternion properties employed in this dissertation is stated in the 
Appendix. Note, that in Equations (104)—(106) z e [0; l] an abstract argument is used 

instead of time. This allows us to exploit certain attributes of the Bezier polynomials and 
define properties at the beginning and endpoints. For example, the analytic expressions 

of the 5th order Bezier Polynomial, 5 (r) > are as follows: 

P l5 (r) = r 5 -5r 4 +10r 3 -10r 2 +5r, 

^ 25 ( r ) = -4r 5 +15r 4 -20r 3 +10r 2 , 

^ 35 (r) = 6r 5 -15r 4 + 10r 3 , (107) 

A, 5 (0 = -4r 5 +5r 4 , 

^ 55 (r) = r 5 . 

In this case, 

q(z-) = q 0 exp 5 (r)) exp (<» 2 /? 2 5 (r)) exp (to 3 ^ 3 5 (r)) exp (tb 4 /? 4 5 (r)) exp (<b 5 /3 55 (r)).(108) 

A plot of the 5th order polynomials is shown in Figure 74. 
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Figure 74. Value of 5th order Bezier Polynomials with respect to the virtual arc. 
These expressions have the favorable properties: 

A 5 (0) = 0 and P i 5 (1) = 1, for / = 1,...,5; (109) 

A 5 (0 ) = As(1 ) = 5, A5(1)=A,(0) = 0, 

A 5 (0) = A 5 0) = 0, for/= 2,3,4; 

A, 5 (0) = A, 5 (l) = 20, A 5 (1) = A. 5 (0) = 0, 

A. 5 (0) = A, 5 (l) = -20, A,(1) = A. 5 (0) = 0, (111) 

A, 5 (0) = A 5 (1) = o, 

which fix the value of the polynomial and its derivatives at the endpoints specified by 
values of re [0; T f ] where we set r f = I that will act as an abstract virtual time argument 

to exploit the above properties in Equations (109)—(111). The values of b) ( are 
coefficients defined by the following relation: 

= Mq^q,), for / = 1,..., 5 . (112) 
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Here, q are the constant column vectors that act as control points and ci, represents a 
constant augmented angular velocity vector based on Equation (112). 

The prevailing idea is that at q(r = 0) = q 0 and q(r = l) = q 5 , where q 0 and q 5 

can be fixed such that q 0 =q(t 0 ) and q 5 = c] (/). This also results in a straightforward 

calculation of higher order derivatives of the quaternion curve with respect to the virtual 
domain argument t . Results for the first derivative were presented in Kim for a 3rd order 
Bezier polynomial (Kim et al. 1995). The results of the 2nd order derivative of a 5th 
order Bezier based quaternion is shown in Equation (113): 
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^3 | ^3 


= q 0 expOfl^ (r ))(«,/?, 5 (r)) 2 ex p(w 2 A, 5 ( r )) exp(d) 3 /? 3 5 (r)) exp(oj 4 /? 4 5 (r))exp(to 5 /? 5 5 (r )) 

r=0 

+ q 0 exp(«! P x 5 (^))(W I ^ ] 5 (r)) exp(« 2 /? 2 5 (r)) exp(<» 3 /? 3 5 (r)) exp (<s 4 /? 4 5 (r)) exp(<» 5 /? 5 5 (r)) 

+ q 0 exp(o) l /i l5 (r))( 03 ] /? l5 (r)) exp(o) 2 /i 25 (t))(6)JJ 25 (r)) exp(e) 3 /i 35 (r)) exp(o) 4 /\-(r)) exp(o') 5 /\ 5 (r)) 

+ q 0 exp(«,A, 5 (r))((aj 15 (r)) exp ((b 7 J 25 (r)) exp(<» 3 /? 3 5 (r))(<» 3 /? 3 5 (r)) exp(<» 4 /? 4 5 (r)) exp(<» 5 /? 5 5 (r)) 

+ q 0 expCra^ 5 (>))(<»,/?! 5 (r)) exp(« 2 /? 2 5 (r)) exp(« 3 /? 3 5 (r)) exp(« 4 /? 4 5 (r))(<» 4 /? 4 5 (r)) exp(« 5 /? 5 5 (r)) 

+ q 0 expCw^^ (r))(w 1 ^ 15 (r)) exp(<» 2 /? 2 5 (r)) exp(w 3 ^ 3 5 (r)) exp(<» 4 /? 4 5 (r)) exp(<» 5 /? 5 5 (r))(<» 5 /? 5 5 (r)) 

+ ^(q° expCw^j 5 (r)) exp(<» 2 /? 2 5 (t))(g> 7 J 25 (r)) exp(<» 3 /? 3 5 (r)) exp(<» 4 /? 4 5 (r)) exp(<» 5 /? 5 5 (r)) J 
+ ^(q° expCWj^ 5 (r)) exp(d> 2 A, 5 0)) exp(<» 3 /? 3 5 (r))(<» 3 /? 3 5 (r)) exp(<» 4 /? 4 5 (r)) exp(<» 5 /? 5 5 (r))j 
+ ^(q° exp(©j # 5 (r)) exp(<» 2 /? 2 5 (r)) exp(<» 3 /? 3 5 (r)) exp(« 4 /? 4 5 (t))(g)J 45 (r)) exp(<» 5 /? 5 5 (r))j 
+ ^ (q° exp(© i; 0 li5 (r)) exp(« 2 /? 2 5 (r)) exp(« 3 /? 3 5 (r)) exp(« 4 /? 4 5 (r)) exp(« 5 /? 5 5 (r))(<» 5 /? 5 5 (r))j. (113) 
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Differentials of quaternion trajectories based on varying order Bezier polynomials 
can be analytically calculated by proper application of the chain rule. The polynomial 
was expanded from a 3rd to a 5th order in order to fix the first and second derivative of 
the quaternion function at the endpoints. By applying Equations (105)—(109), therefore 

exploiting the property that certain terms /?. 5 , /?. 5 and /) 5 are equal to zero, derivative 
values at the endpoints can be calculated by the simple expression: 


dq 

dr r=0 


~ dq 
5q 0 «i, — 
dr 


= 5q 5 w 5 


T =1 


(114) 


Note, in Equation (114), the first derivative, how the parameter tb ( is related to the 
angular velocity at the endpoints. In similar fashion, 


d 2 q 

~d? 


Ir=0 


= ^Oq,,®! + 25q 0 d) 1 2 + 20q 0 d) 2 , 


d 2 q 

dr 2 


It=1 


= -20q 4 w 4 q 4 'q 5 +20q 5 d> 5 +25q 5 ra 5 2 . 


(115) 


2. Mapping from the Virtual Arc to the Time Domain 

Now that the trajectory is set using a virtual domain, a mapping must be 
employed to convert this trajectory into a time dependent one. To do this, a speed 
factor, A , is defined that maps the points on the trajectory from the virtual domain to the 
time domain, therefore defining the final time of the maneuver: 


A = 


dr 
~dt ’ 


dr 

T 


(116) 


The Speed factor can be a constant parameter that simply stretches or shrinks time 
evenly or it can take the shape of more complex function. For this application, A(r) is 
restricted to a function that contains a reduced number of varied parameters, developed 
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for this dissertation, that still allow for the speeding up and slowing down along a 
trajectory. The representation is shown in Equation (117). 

A(t) = A 0 + At 2 + (1 - r) 2 f? + (l - (1 - r) 2 ) C + (1 - r 2 )D (117) 

Obviously, more complex structures of A(z) will provide more flexibility in the 
trajectory. One main idea of Equation (117) is to keep A(z) positive for all time, as 
A(t )<0 would imply time marching backwards and a(t) =0 signifies that time has 
stopped. Also, keeping the speed factor of a form such that an analytic integral to 
Equation (117) exists not only provides computational efficiency and an accurate 
integration to the minimum time PI but also provides a continuous mapping from the 
virtual domain to the time domain. Alternatively stated, a continuous control history is 
available, whose resolution does not suffer from a limited number of node points. This 
attribute is later revisited in Chapter VI.D.2. 

Although a speed vector of the form Equation (117) does not allow matching the 
optimal minimum time solutions exactly, varying the parameters contained within A(z) 
(Ag, A, B, C, and D) still allows variation of the speed along the trajectory defined by 
Equation (104) to produce feasible and easy to track suboptimal solutions. 

3. Inverting the Dynamics 

As a result of the mapping from virtual to time domain, the expression for the 
differential of q with respect to time is: 

4 = = = (ii8) 

dt dr dt 2 dr dz dz 

Now, if the trajectory in the virtual domain, q(r), is specified, along with the 
speed trajectory, A(z), the resulting trajectory of ,q(/), as well as its higher order 
derivatives, can be analytically expressed and mapped to the time domain. 
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Inverting kinematic equations (100) and differentiating the result yields analytical 
expressions for the angular velocity and angular acceleration: 

©(0 = 2q _1 (0q(0» (119) 

a(0 = q _1 (02q(0-q(0<o(0- 

The torque history needed to follow such a trajectory is calculated by inverting Equation 
(99): 

T x ( t ) = a x (0 + (/ 33 - 1 22 )co y (; t)co z (t), 

T y (t) = a y (t) + (I u -I 33 )oj t (t)co z (t), (120) 

T z ( t ) = a z (0 + (/ 22 - /,, )cu v (0®, (0- 


4. Matching Endpoint Conditions 

From the preceding equations, a quaternion history can be developed based on the 
Bezier polynomial that satisfies predefined beginning and ending quaternion values as 
well as setting the angular velocity and angular acceleration at the endpoints. The desired 
angular velocity and acceleration dictate the quaternion derivatives at the endpoints to be 
as follows: 


q(V) = q(V)t°(h.)’ 

q(tf) = q(t f )o)(t f ), 

q(* 0 )=q( ? o) 2a (^o)+q(^o) o >(^oX 

q (tf) = q(t f )2a(t f ) + q (t f )(o(t f ). 


( 121 ) 


Based on the properties of the 5th order Bezier polynomial, the coefficients of the 
quaternion expression are then calculated by: 


=^(t 0 ), (122) 

w 5 =^ t0 (^/) 5 ( 123 ) 


to 2 


q 0 1 q(to) + 20c)i 
20 


-25«! 2 


(124) 
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(125) 


to 4 


q 4 -1 (q(^/) + 2 °qs t Qi - 25q 5 co 1 2 )q.-’q, 
-20 


Here, q is computed using the second equation in Equation (119) and the 
complementary q ( parameters defined as: 


and 


Finally, 


qo=qOo) = q(OLo’ 

q 5 =q(^ / ) = q(«’)U» 

qi = q 0 e x P(tOi), 

q 4 = q 5 exp(d) 5 ) -1 . 


q 2 = q! exp(tb 2 ), 

q 3 = q 5 (exp(d> 4 )exp(d> 5 )) 


® 3 = ln (q2~'q 3 )- 


(126) 


(127) 


(128) 


The resulting benefit of this laborious fonnulation is that the attitude trajectory 
history of a 4x1 q vector, that satisfies the constraints of a unit quaternion, can be 
specified by a reduced set of parameters. These parameters are the initial and final 
conditions on the quaternion itself, as well as values of angular velocity and angular 
acceleration at those endpoints. 


5. Increasing the Polynomial Order 

More flexibility in the trajectory can be obtained by increasing the order of the 
Bezier polynomial used in the basis function. For the case of a 7th order polynomial, the 
same structure as Equation (104) is employed but now with n=l. The resulting 7th order 
polynomials are shown in Figure 75. 
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Figure 75. Value of 7th order Bezier Polynomials with respect to the virtual arc. 


This leads to the introduction of q 6 ,q 7 , d> 6 and d> 7 , which are defined to be 

consistent with previous definitions from Equations (104) and (112). If the values of 
orientation and angular velocity at the endpoints are set, endpoint conditions of angular 
jerk as well as angular acceleration can be used as varied parameters. Setting a specified 
(low) value for the initial and final jerk can be critical for slewing maneuvers of flexible 
spacecraft, in order to avoid excitation of structural modes. To do this, the third 
derivative of the output vector is calculated as: 


<i 3 q , 2 <i 2 q , dA , d 2 q dA 2 dx\ d 1 /i . 2 dqfdA 

- A H-— 2 A A H--- A H- —A +• 


dr 


dr dr dr dr 


dr dr" 


dr 


dr 


A. 


(129) 


The new expressions for the virtual derivatives are also recalculated, which are analogous 
to Equations (114) and (115) but with the addition of a third derivative to accommodate 
the change between 5th and 7th order polynomial: 
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d 3 q 


dr 3 


= 343q 0 ®i -882q 0 tof -420q 0 to 2 +210q 0 d) 1 +210q 0 d) 3 +882q 0 to 1 d) 2 , 


r=0 


d 3 q 


dr 


( 130 ) 


= 343q 7 d) 7 -882q 6 d) 6 e) 7 q 6 1 q 7 + 210to 5 q 7 +210q 7 to 7 +882q 7 d) 7 


T=\ 


+ 882q 0 to 1 d) 2 -420to 6 q 7 ; 


d 2 q 


dr 
d 2 q 


= 49q 0 d) 1 2 + 42q 0 d) 2 - 42q 0 d> 1 , 


T= 0 


dr 


= 42q 7 d) 7 + 49q 7 d) 7 - 42co 6 q 7 ; 


r=l 


(131) 


iiq 

dr r=0 

iiq 

dr r=1 


= 7q 0 ©i, 

= 7q 7 w 7 . 


(132) 


New values for the constants that fix the initial conditions of the quaternion 
trajectory can be calculated similar to the 5th order polynomial, except and extra step 
needs to be taken to accommodate the third derivative of q . 

D. SOLVING THE PROBLEM USING IDVD METHOD 

This section presents the results of using IDVD method with two different 
parameterizations to obtain the minimum time solutions of the problem posed in Chapter 
III.A. The Matlab function fmincon was used to optimize the trajectory while varying 
either the angular acceleration (for the 5th order polynomial) or the angular acceleration 
and jerk (for the 7th order polynomial) at the endpoints of the trajectory. Together with 
five varied parameters for X(t), Equation (117), it yields 11 varied parameters for the 5th 
order polynomial approximation and 17—for the 7th order polynomial approximation. 
The constraints are that the resulting control must obey Equation (102) and T(t) cannot 
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be < 0 at any instant. Initial guesses for the angular acceleration and jerk are equal to 
zero and the initial guess of 10‘ 4 is used for all parameters contained in A (t ). 

1. IDVD Solutions Varying 2nd Derivative Results 

Figures 76-79 present the results obtained when applying the IDVD with a 
quaternion based on a 5th order Bezier polynomial (compare it with the GPOPS solution 
presented in Figures 66-68), allowing variation of the 2nd derivative at the endpoints. 
The solution was evaluated using 100 nodes (although, as opposed to GPOPS because the 
solution is analytic, it would be no difference running it for a larger or smaller number of 
nodes), resulted in slightly larger t , but took significantly less time to compute at only 

10.0 seconds. 


-1-1-1-1-1- 

♦ ii 

-q 2 
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llqll 

. 

l_l_l_1_1_l 
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Figure 76. Case 1 (IDVD 5th order): time history of the quaternion. 
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Figure 77. Case 1 (IDVD 5th order): time histories of the angular velocity and 

accelerations. 



Figure 78. Case 1 (IDVD 5th order): time history of controlling torques. 

The major difference compared to the GPOPS solution is that the controls do not 
have a bang-bang nature anymore. Again, it was done intentionally by the choice of the 
quaternion parameterization. It occurs that when implemented in the real time controller, 
these controls will be easier to track. Also, having different controlling torques at the 
endpoints means having different angular accelerations. While for the GPOPS solution, 
the terminal angular accelerations are at the mercy of the optimization routine using 
IDVD allows matching them with the current accelerations, which also makes the control 
algorithm more robust. Figure 78 shows the speed factor, the key element in matching 
the virtual and time domains. 
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Figure 79. Case 1 (IDVD 5th order): mapping the virtual and time. 

Figure 80 shows an outline of the slew in inertial space, clearly showing that it is not an 
eigenaxis maneuver. 



Figure 80. 


Case 1 (IDVD 5th order): principal axis outline of 180° slew. 
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The resulting solution for the same scenario using 25 nodes is shown in Figures 81-84. 


Figure 81. 


Figure 82. 


Figure 83. 



0 0.5 1 1.5 2 2.5 3 3.5 


Time, s 

Case 1 (IDVD 5th order): time history of the quaternion using 25 nodes. 




Case 1 (IDVD 5th order): time histories of the angular velocity and acceleration 

using 25 nodes. 



Case 1 (IDVD 5th order): time history of controlling torques using 25 nodes. 
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Figure 84. Case 1 (IDVD 5th order): mapping the virtual and time domains using 25 

nodes. 

The rate profiles for the solution of 25, 50, and 100 nodes are shown in Figures 85 
and 86. Each case results in a smooth control solution regardless of the number of nodes 
chosen. This is due to the construction of the quaternion history as well as the controls. 
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Figure 85. 


Figure 86. 




Case 1 (IDVD 5th order): comparison of angular velocity time histories, 
obtained for 25 nodes (top) and 50 nodes (bottom). 




Case 1 (IDVD 5th order): comparison of torque time histories, obtained for 25 
nodes (top) and 50 nodes (middle). 
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As in the case of GPOPS solution for Case 2, the nonsymmetrical inertia matrix 
causes a certain changes as compared to the symmetric matrix solution. The 100-node 
IDVD solution in this case results in 4.767s maneuver and requires about a minute to 
compute (see Figures 87-90). 



Figure 87. Case 2 (IDVD 5th order): Time history of the quaternion. 



0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 



Figure 88. Case 2 (IDVD 5th order): time histories of the angular velocity and 

acceleration. 
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Figure 89. Case 2 (IDVD 5th order): time history of controlling torques. 



Figure 90. Case 2 (IDVD 5th order): mapping the virtual and time domains. 

2. IDVD Solutions Varying 2nd-3rd Derivative Results 

For the sake of comparison, Figures 91-94 present the solution of the same 
problem using a quaternion based on a 7th order Bezier polynomial, allowing variation of 
both the 3rd and 4th derivative at the endpoints. Although this slightly improves the PI, it 
drastically increases the computational time. The following section addresses this issue in 
more detail. 
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Figure 91. 




Case 1 and 2 (IDVD 7th order): time history of the quaternions Case 1 (top) 

and Case 2 (bottom). 
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Figure 92. 
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Case 1 and 2 (IDVD 7th order): time hi s tories of the angular velocity and 
acceleration for Case 1 (top) and Case 2 (bottom). 
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Figure 93. 




Case 1 and 2 (IDVD 7th order): time history of controlling torques for Case 1 

(top) and Case 2 (bottom). 
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Figure 94. Case 1 and 2 (IDVD 7th order): mapping the virtual and time domains for Case 

1 (top) and Case 2 (bottom) solutions. 

E. RESULTS AND COMPARISONS 

This section presents a comparison of the results obtained using the IDVD method 
with those of the GPOPS method. It disregards the fact that the results obtained with 
GPOPS for low number of nodes are infeasible but rather concentrates on the 
computational advantages the IDVD approach has for any number of intennediate points 
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(nodes in the case of GPOPS). To start with, Tables 17 and 18 summarize the 180° rest- 
to-rest slew maneuver solutions for symmetric and asymmetric inertia matrix obtained 
using GPOPS and IDVD, as discussed in Chapters IV.C and IV.D. 

In these tables, all results are compared against the readily available eigenaxis 
maneuver solution. First, it is shown that the true optimal solution, obtained offline 
(Bilimoria and Wie 1993), which is not an eigenaxis rotation, provides about 8.5% and 
11% improvement of the PI, time of maneuver, for Case 1 and Case 2, respectively. 


As seen from Tables 17 and 18, the GPOPS solution converging to one of the 
equally optimal solutions (if at all), assures about the same gain in the PI as the truly 
optimal one. But, again, it takes significant computational time. Specifically, the 100- 
node solution that does converge and assures a smooth controls history takes about an 
hour to converge. 


Trajectory Generation 
Method 

Nodes 

Computational 
Time (sec) 

cost (ff) 

% 

improvement 

Eigenaxis 

N/A 

N/A 

3.5449 

O 

'sP 

Optimal (BilimoriaAA/ie) 

N/A 

N/A 

3.2431 

8.51 % 

GPOPS 

25 

15.5 

3.2573 

8.11% 

GPOPS 

50 

200.5 

3.2859 

7.31 % 

GPOPS 

100 

5962.8 

3.2430 

8.52% 

IDVD 5th order 

25 

3.7 

3.4289 

3.27% 

IDVD 5th order 

50 

4.8 

3.4373 

3.04% 

IDVD 5th order 

100 

10.0 

3.4382 

3.01 % 

IDVD 7th order 

25 

44.0 

3.3756 

4.78% 

IDVD 7th order 

50 

69.8 

3.3769 

4.74% 

IDVD 7th order 

100 

121.2 

3.3776 

4.72% 


Table 17. The 180° rest-to-rest slew maneuver about the z-axis, symmetric inertia (Case 1). 
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Trajectory Generation 
Method 

Nodes 

Computational 
Time (sec) 

cost (tf) 

% 

improvement 

Eigenaxis 

N/A 

~0 

5.0133 

~0% 

GPOPS 

25 

22.2 

4.4609 

11.02% 

GPOPS 

50 

520.1 

4.4499 

11.24% 

GPOPS 

100 

1893.9 

4.4528 

11.18% 

IDVD 5th order 

25 

4.8 

4.7857 

4.54% 

IDVD 5th order 

50 

7.6 

4.7861 

4.53% 

IDVD 5th order 

100 

9.5 

4.7864 

4.53% 

IDVD 7th order 

25 

24.3 

4.6768 

6.71 % 

IDVD 7th order 

50 

41.2 

4.6846 

6.56% 

IDVD 7th order 

100 

80.4 

4.6869 

6.51 % 


Table 18. The 180° rest-to-rest slew maneuver about the z-axis, asymmetric inertia (Case 2). 


As discussed in the previous section, the IDVD solution has a much more robust 
performance, allowing computing the same type of maneuvers just in a few seconds as 
opposed to hours. If an executable optimization library is implemented, the IDVD 
method produces solutions in fractions of a second (Yakimenko and Siegers 2009). Of 
course, some of the optimality (PI value) is sacrificed. On the positive side, the solution 
is always feasible and smooth for any number of computational points, and can be 
brought closer to the GPOPS solution (in terms of the value of a PI) by increasing the 
number of varied parameters (the order of the quaternion approximation polynomial). 
Furthermore, since the resulting IDVD solution is analytic in nature, increasing the node 
points, after a solution is obtained, is a trivial evaluation of an analytic expression. 

As discussed in Chapter IV.B, this solution features a bang-bang control, i.e., does 
not account for controllers’ dynamics, and therefore can still not be used onboard as is. 
On the other hand, the always-feasible and ready-to-go IDVD solution (employing as low 
as say 25 computational points) can be produced much faster, but surrenders up to 2 A of 
its gain as compared to that of the GPOPS solution (about Vi for the 7th order 
approximation). 

Tables 19 and 20 present similar data for the 90° rest-to-rest slew maneuver. 
While GPOPS provides about 3% gain compared to a simple eigenaxis slew solution, the 
IDVD method has almost no advantage or may be even worse if using a 5th order 

quaternion approximation. This is because, at some point, the gains made by having the 
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ability to calculate a bang-bang solution (with discontinuous controls) outweighs the 
benefits of having the ability to exploit the gains made by the effect of slewing off the 
eigenaxis. All major conclusions, however, remain the same. 


Trajectory Generation 
Method 

Nodes 

Computational 
Time (sec) 

cost (tf) 

% 

improvement 

Eigenaxis 

N/A 

N/A 

2.5066 

~0% 

GPOPS 

25 

20.7 

2.4336 

2.91 % 

GPOPS 

50 

120.0 

2.4332 

2.93% 

GPOPS 

100 

236.2* 

2.4296 

3.07% 

IDVD 5th order 

25 

5.8 

2.5654 

-2.34% 

IDVD 5th order 

50 

7.2 

2.5666 

-2.39% 

IDVD 5th order 

100 

9.4 

2.5671 

-2.41 % 

IDVD 7th order 

25 

39.7 

2.5043 

0.09% 

IDVD 7th order 

50 

58.9 

2.5058 

0.03% 

IDVD 7th order 

100 

77.3 

2.5058 

0.03% 


Table 19. The 90° rest-to-rest slew maneuver about the z-axis, symmetric inertia (Case 1). 


Tables 21 and 22 compare GPOPS and IDVD solutions for such case, when other 
sets of nonzero boundary conditions were explored as well, and proved to maintain the 
same trends. In this table, all results are compared against a valid 100-node GPOPS 
solution. As seen, the GPOPS solutions with a lesser number of nodes produce 
somewhat infeasible solutions, meaning that they cannot be implemented in the control 
scheme explicitly. The IDVD solutions may yield to GPOPS as much as about 4% with 
respect to the PI, but again are produced much faster. Furthermore, as shown in Figures 
95 and 96, the 90° and 180 0 rest-to-rest slew maneuvers with zero boundary rates feature 
multiple equally optimal solutions, so that both GPOPS and IDVD solutions converge to 
different solutions when changing the number of nodes (GPOPS) / computational points 
(IDVD). In contrast, for the case of nonzero boundary rates, they all converge to the 
same solution,as shown in Figure 97. 
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Trajectory Generation 
Method 

Nodes 

Computational 
Time (sec) 

cost (ff) 

% 

improvement 

Eigenaxis 

N/A 

~0 

3.5449 

~0% 

GPOPS 

25 

17.4 

3.4450 

2.82% 

GPOPS 

50 

168.0 

3.4408 

2.94% 

GPOPS 

100 

1432.1 

3.4430 

2.87% 

IDVD 5th order 

25 

5.9 

3.6277 

-2.33% 

IDVD 5th order 

50 

9.8 

3.6307 

-2.42% 

IDVD 5th order 

100 

17.3 

3.6316 

-2.44% 

IDVD 7th order 

25 

41.5 

3.5373 

0.21 % 

IDVD 7th order 

50 

80.4 

3.5389 

0.17% 

IDVD 7th order 

100 

158.8 

3.5394 

0.16% 


Table 20. The 90° rest-to-rest slew maneuver about the z-axis, asymmetric inertia (Case 2). 


It should be noted that in practice, the direct methods would likely be used in 
situations where the end-conditions (angular rates, accelerations) of the slew are specified 
and not equal to zero (to meet mission requirements of matching attitude rates of a 
tumbling vehicle, for example). For this case, no simple eigenaxis slew solution exists 
and therefore any solution produced on-line would be good. The next chapter of this 
dissertation exploits the IDVD and coupling with translational maneuvers. 


Trajectory Generation 
Method 

Nodes 

Computational 
Time (sec) 

cost (fi) 

% within 
GPOPS 100 
nodes 

GPOPS 

25 

48.3 

2.4028 

0.07% 

GPOPS 

50 

333.9 

2.4016 

0.02% 

GPOPS 

100 

1182.6 

2.4011 

0.00% 

IDVD 5th order 

25 

5.9 

2.5437 

5.94% 

IDVD 5th order 

50 

5.6 

2.5450 

5.99% 

IDVD 5th order 

100 

6.1 

2.5451 

6.00% 

IDVD 7th order 

25 

27.9 

2.4885 

3.64% 

IDVD 7th order 

50 

41.5 

2.4895 

3.68% 

IDVD 7th order 

100 

90.6 

2.4896 

3.69% 


Table 21. The 90° maneuver for symmetric inertia (Case 1) and nonzero boundary rates, 

e> 0 = ~(o f = [0.1 0.1 0.l] T . 
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Trajectory Generation 
Method 

Nodes 

Computational 
Time (sec) 

cost (t^) 

% within 
GPOPS 100 
nodes 

GPOPS 

25 

11.7 

2.9054 

0.03% 

GPOPS 

50 

111.3 

2.9047 

0.01 % 

GPOPS* 

100 

164.3 

2.9045 

0.00% 

IDVD 5th order 

25 

6.3 

3.0060 

3.49% 

IDVD 5th order 

50 

5.5 

3.0064 

3.51 % 

IDVD 5th order 

100 

8.3 

3.0066 

3.52% 

IDVD 7th order 

25 

31.4 

2.9559 

1.77% 

IDVD 7th order 

50 

60.1 

2.9568 

1.80% 

IDVD 7th order 

100 

106.5 

2.9571 

1.81 % 


Table 22. 


The 90° maneuver for symmetric inertia (Case 1) and nonzero boundary rates, 

CO 0 = ~(O f = — [l 1 l] T . 
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Figure 95. Case 1: 180° rest-to-rest slew profile as projected onto the x-y inertial plane. 
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Figure 96. Case 1: 90° rest-to-rest slew profile as projected onto the x-y inertial plane. 
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Figure 97. Case 1: 90° slew profde, co 0 = -to f = [0.1 0.1 0. l] T , as projected onto the x-y inertial plane. 
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V. RAPID ATTITUDE AND TRANSLATIONAL TRAJECTORY 

GENERATION 


Inverse Dynamics in the Virtual Domain exploits the ability to take a polynomial 
of sufficiently high order and match desired values, and potentially derivative values, at 
its endpoints to generate a trajectory. This is done by representing each position variable 
of the trajectory of the close approach maneuver in the orbital frame by a summation of 
polynomials and trigonometric basis functions. While some endpoint conditions are 
specified based on the problem statement, higher order derivative of the endpoints can 
function as parameters that may be varied to optimize the trajectory based on a chosen PI. 
The remainder of this chapter explains the process in detail and applies the IDVD method 
to the benchmark problem stated in Chapter III. 

A. INVERSE DYNAMICS IN THE VIRTUAL DOMAIN FORMULATION 

First, a basis function for the translational trajectory is defined by Equation (133). 


r 7i 


n 


p( r) = a 0 +a l r + a 2 T^ + a 2 r +b 0 sin(—r) + c 0 cos(—r ) + b l sin(^r) + Cj cos(^r) . (133) 


This leads to 3 separate equations (and separate set of coefficients) for the x, y, and z 
position of the chaser vehicle each based on the form of the basis function described by 
Equation (133). It should be reinforced that the trajectory described by Equation (133) 
with higher order derivatives shown by Equations (134)—(135), are constructed with 

respect to a virtual argument, r e ^(kz^J, and not time. This specifies the trajectory in 

the spatial domain, but lets the speed the trajectory is traversed to be varied as well. 


P = 


dp 

dv 


f 7T 77 77 77 ^ 

<7j + la 2 r + 3a 3 r +b 0 —cos (—r)-c 0 —sm(—T) 


+b { 7i cos(^t) - cpi sin(^r) 


J 


dp 
~d~T Y 


^ TC > . 7t Tt 1 77 ^ 

a 2 +6a 3 r-h 0 — sin(— r) + c 0 — cos(— r) 
.2 


-b x 7T sin(^r) -cpt cos {nx') 


J 


(134) 


( 135 ) 
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The plots of the basis functions considered are shown in Figure 98: 
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Figure 98. Basis functions considered for translational trajectory generation. 


1. Mapping the Virtual Arc and Specifying Endpoint Values 

The same structure of the virtual time argument is used as in Chapter IV.B.2, 
Equation (117). This allows for 5 varied parameters that determine the history of the 

speed factor, ^, over the duration of the maneuver. The endpoint values of the basis 
functions are simply set to the beginning and ending values of the maneuver in the spatial 
domain. The derivative values at the endpoints and their derivatives in the virtual domain 
are calculated based on the time derivatives at the endpoints and speed factor shown in 
Equations (136) and (137). 


dp jit f ) 
dx 


= Pi(tf)/A(t f ), 


dPi(Jo) 

dx 


for i = x,y,z 


Pi(t 0 )/A(t 0 ) 


(136) 
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and 


d 2 Pi(t f ) dA(t f ) 

d PM o) \ ■ u \ dMt 0 ) 




(137) 


dr 2 


■ = (Pi(t 0 )-Pi(t 0 ) 


dt 


)U(t 0 y- 


In order to specify the 2nd derivatives at the endpoints, a basis function with a reduction 
in the number of coefficients is used (Yakimenko 2000; Yakimenko and Siegers 2009): 


,7t 


p(r) = a 0 + a x z + apr + c 0 cos(— z ) + b \ sin(^r) + c l cos(ttt) . 


(138) 


The resulting coefficients are then be derived according to the algebraic relation in 
Equation (139). 
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For higher order basis functions that are capable of specifying jerk at the 
endpoints, the basis function is specified as Equation (133) and the coefficients are 
obtained by solving the algebraic relation in Equation (140). 
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2. Implementation of Quaternion Based Attitude Representation 

The attitude implementation of IDVD for the close rendezvous problem described 
in Chapter III is done similar to Chapter IV.B. The proper Bezier polynomial is chosen 
such that the desired boundary conditions can be met as well as allowing variation of the 
next highest derivative. One major difference is that the attitude trajectory is formulated 
in the orbital frame instead of the inertial. The angular rates with respect to the inertial 
frame are then computed, based on Equation (141), and used to detennine the resulting 
control torques in the body frame. 

(D x — co x + 2(q l q i + q 2 q 4 )fi 

a>y = °G>y + 2 (<h^ ~ c h c U) D - ( 141 ) 

co, — co, + 2(q 4 — q l — q 2 + q 2 )Q 

3. Alternative Euler Angle Formulation for Attitude Representation 

Considering an Euler angle representation of a 2-1-3 with respect to the orbital 

reference frame rotation can be expressed by the angular displacements an< ^ 

pertaining to angular displacements about the y, x and then z body axis respectively (Sidi 
1997). The angular rates with respect to the orbit frame can then be expressed as: 

°co x = (j) cos {y/) + 6 sin(^) cos(^) 

° co y = -(j) sin(^) + 6 cos(^) cos(^) (142) 

°co z =\j/ -0 sin(^) 

The angular rates with respect to the inertial frame can then be expressed as: 

co x =0cos(^) + ^sin(^/')cos(0) + 2 (< 7 1 <73 +q 2 q 4 )Cl 

co y = -^sin(^-) + ^cos(^)cos(^) + 2(^ 2 ^3 -q x q A )Cl (143) 

co, = \j/ - 6 sin(^) + 2{q 4 -q , 2 -q 2 2 +q 2 )Cl 

The resulting angular accelerations, with respect to the inertial frame expressed in the 
body frame becomes: 
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a x = (f) cos (y/) ~ (frys sin (y/) + 6 sin(^) cos(^) + 6\j/ cos {(/)) cos(^) 

- 6<j) sin ((f)) sin(^) + 2(q 3 q x + q x q 3 + q A q 2 + q 2 q A )a> 0 

a =-^sin(^-)-^cos(^) + ^cos(^)cos(^)-^cos(^)sin(^) (144) 

- Q(j) sin(^) cos(^) + 2(q 2 q 3 + q 3 q 2 - q 4 q x - q x q 4 )co o 
a z =\j/- 6(j) cos {f) - 6 sin(^) + 2 (q 4 q 4 -q x q x - q 2 q 2 + q 2 q 2 )a> 0 

This equates to the required inertial torque expressed in the body frame dictated by: 

T x =I n a x +{! i3 -I 72 )a> z <D y 

Ty = / 22« +(/|, - 4 K®, ( 145 ) 

T z — Ar^i ^ (^22 i )^,(O x 

This formulation allows complete description of the attitude trajectory in the orbital and 
inertial frames and the associated rates and inertial torques by specifying initial and final 
conditions on 9, t//, and (f) as well as their higher order derivatives. 


B. SOLVING THE IDVD PROBLEM 

For the close rendezvous problem with the IDVD formulation, the Sparse 
NOlinear OPTimizer (SNOPT) (Gill, Murray, and Saunders 1996) was used to solve for 
the resulting trajectory. For this implementation, the objective function was the PI based 
on the specific type of cost function. The problem constraints were programmed as: 


control constraint: -u . < u < u mav 

min nidA 

path constraint: (x 2 + y 2 +z 2 )-r 2 > 0 

time constraints: A(r) > 0 V r 


(146) 


with the limits on u and value of r set to those used in Chapter III. It is reinforced that, 
with the IDVD method, the dynamic constraints (equations of motion) of the problem and 
desired endpoint conditions are always satisfied exactly at every iteration. If the problem 
is not stated such that a solution does not exist, judicious (or sometimes common sense) 
selection of the initial guess on the varied parameters, conservative guesses on the 
endpoint derivatives and speed factor coefficients, will lead to always having a feasible 
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solution to implement. This is not the case with pseudospectral methods, since at any 
given iteration, the kinematic and dynamic constraints in the equations of motion are not 
guaranteed to be satisfied and may result in an infeasible solution. 

C. COMPUTATIONAL RESULTS 

Both IDVD methods of matching endpoint translational and angular position and 
velocities are used to generate solutions. The 1st set of endpoint parameters is considered 
to be the position and the 2nd set of endpoint parameters is the velocity. The IDVD (3rd) 
method refers to varying the set of 3rd parameters, accelerations, at the endpoints while 
the IDVD (3rd-4th) method refers to varying both the third and fourth parameters, 
accelerations and jerk. 

1. Minimum Energy Cost Varying 3rd Parameter Set Results 

The results for using the IDVD (3rd) method for a minimum energy cost, 
Equation (15), are shown next. They are based on the translational basis Equation (139) 
and the quaternion attitude formulation employing the 5th order Bezier polynomial. The 
result is suboptimal when compared to the infinite dimensional optimal control problem 
solution, but is optimal based on the additional constraint of only using the specified 
polynomial set of basis functions. This reduces the number of varied parameters to 17 
(the third derivative of each state plus the coefficients of the speed factor), resulting in a 
computational time of 12.2818 seconds. The final PI based on the minimum quadratic- 
control cost is J = 0.2461. Figure 99 shows the 3D representation of the trajectory in 
Figure 100 shows the planar projection of the same trajectory. The control solution is 
shown in Figure 101, which is smooth over the interval, consistent with the optimal 
solution found in Chapter III.C.2. Figures 102 and 103 show the time history of the 
endpoint discrepancies and verify that the endpoint conditions were met at the final time. 
The solution using this specific formulation and cost is highly attractive due to it’s 
inherent low fuel cost and the attribute that the solution does not demand that the 
actuators be saturated at any given instant, unlike the bang-bang and bang-off-bang 
nature of minimum time and minimum fuel cost functions respectively. Specifically, this 
noncontrol-saturated, smooth trajectory is desirable as the spacecraft is in the final stages 
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of matching docking points. The speed factor, a representative measure of how fast the 
chaser is moving along the spatial trajectory, is shown in Figure 104. 



Figure 99. 


Minimum energy solution (IDVD 3rd): 3D optimal rendezvous trajectory. 
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Figure 100. Minimum energy solution (IDVD 3rd): 2D planar projection of the trajectory. 
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Figure 101. Minimum energy solution (IDVD 3rd): control history. 



Figure 102. Minimum energy solution (IDVD 3rd): history of the translational endpoint 

discrepancies. 
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Figure 103. Minimum energy solution (IDVD 3rd): time history of the rotational endpoint 

discrepancies. 
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Figure 104. Minimum energy solution (IDVD 3rd): time history of the speed factor. 
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2. Minimum Time Cost Varying 3rd Parameter Set Results 

The final PI based on the minimum time cost is J = 3.734 seconds. Figure 105 
shows the 3D representation of the trajectory in Figure 106 shows the planar projection of 
the same trajectory. The control solution is shown in Figure 107, which although is 
smooth by construction, approaches the bang-bang nature of an optimal solution 
computed with respect to a minimum time cost function presented in Chapter III.C.l. 
Figures 108 and 109 show the time history of the endpoint discrepancies and verify that 
the endpoint conditions were met at the final time. Figure 110 shows the speed factor, 
which revels that the chaser is speeding up along the spatial trajectory and then slowing 
down as it approaches the final desired states. This type of qualitative behavior is 
inherent to bang-bang, minimum time solutions to general maneuvers which apply max 
control at the beginning to increase speed of the maneuver, then max control at the end to 
decrease speed and arrive at the desired final states. 



Figure 105. Minimum time solution (IDVD 3rd): the 3D optimal rendezvous trajectory. 
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Figure 106. 



Minimum time solution (IDYD 3rd): 2D planar projection of the trajectory. 
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Minimum time solution (IDVD 3rd): control history. 
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Figure 107. 
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Figure 108. 


Minimum time solution (IDVD 3rd): history of the translational endpoint 

discrepancies. 
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Figure 109. Minimum time solution (IDVD 3rd): history of the rotational endpoint 

discrepancies. 
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Figure 110. Minimum time solution (IDVD 3rd): time history of the speed factor. 


3. Minimum Fuel Cost Varying 3rd Parameter Set Results 

The final PI based on the minimum fuel cost is J = 2.8312. Figure 111 shows the 
3D representation of the trajectory in Figure 112 shows the planar projection of the same 
trajectory. The control solution is shown in Figure 113. Upon close inspection, the 
control solution for each component has more action at the beginning and ends of the 
maneuver, exhibiting qualities of a bang-off-bang structure consistent with the cost 
function presented in Chapter III.C.l. Figures 114 and 115 show the time history of the 
endpoint discrepancies and verify that the endpoint conditions were met at the final time. 
Figure 116 shows the speed factor, which is similar to the minimum energy solution 
presented in the previous section. 
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Figure 112. Minimum fuel solution (IDVD 3rd): 2D planar projection of the trajectory. 
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Figure 113. Minimum fuel solution (IDVD 3rd): control history. 



Figure 114. Minimum fuel solution (IDVD 3rd): history of the translational endpoint 

discrepancies. 
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Figure 115. Minimum fuel solution (IDVD 3rd): history of the rotational endpoint 

discrepancies. 




0.17 

0.16 

0.15 


■g 0.14 

CD 

w 0.13 
0.12 


0 




4 

Time, s 


Figure 116. 


Minimum time solution (IDVD 3rd): time history of the speed factor. 
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4. Minimum Energy Varying 3rd-4th Parameter Set Results 

The problem stated in Chapter V.B.l is also solved using the IDVD (3rd-4th) 
method based on the minimum energy cost function from Equation (15). The results are 
extremely similar to the solution, as shown in the 3-D trajectory of Figure 117, to the 
same problem formulation based on the IDVD (3rd) method, except the increased 
flexibility in the polynomial due to the extra basis functions allows for smoother controls 
at the endpoints, as shown in Figure 118 and later in Figure 123, when compared with 
IDVD (3rd). The solution took 623 seconds to compute and had a PI of J= 0.2458, only a 
0.12% reduction in PI. It should be noted that using this formulation, you can specify 
constraints on the jerk profiles while still using the acceleration profile as your control 
vector. This attribute is based on the polynomial formulation of the trajectory and cannot 
be implemented using pseudospectral techniques. To set constraints on jerk using 
pseudospectral techniques, the control vector would have to be based on jerk and not 
acceleration. 
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Figure 118. Minimum energy solution (IDVD 3rd-4th): control history. 

5. Minimum Time Varying 3rd-4th Parameter Set Results 

The minimum time solution using the IDVD (3rd-4th) method is shown in Figure 
119. It has a PI of J = .3.6582 and took 86.6 seconds to compute. The 3D trajectory 
looks similar to that obtained by IDVD (3rd), but the extra flexibility in the polynomial 
provided by an increased set of basis function allows for shaper transitions in the 
controls, shown in Figure 120, mimicking that of a bang-bang controller. Although the 
computational cost rose to 86.6 seconds, it is still well below the 8+ hour computational 
cost of GPOPS for a 200 node representation. 
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Figure 119. Minimum time solution (IDVD 3rd-4th): 3D optimal rendezvous trajectory. 
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Figure 120. Minimum time solution (IDVD 3rd-4th): control history. 
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6. Minimum Fuel Varying 3rd-4th Parameter Set Results 

The minimum fuel solution obtained using the IDVD (3rd-4th) method is shown 
in Figure 121. Again, the trajectory is very similar to the solution calculated from the 
IDVD (3rd) except the increased flexibility allows for shaper transitions in control at the 
endpoints. The bang-off-bang type of control is evident from examining the control 
histories in Figure 122. The resulting PI is J = 2.7552, taking 405.5 seconds to compute. 
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Figure 122. Minimum fuel solution (IDVD 3rd-4th): control history. 

7. Propagated IDVD Solutions 

As with the solutions obtained by the pseudospectral method in Chapters II and 
III, the solutions for the controls obtained with the IDVD method were propagated using 
the same mechanism as the pseudospectral method stated in Chapter III.C.4. The 
resulting discrepancies in the endpoints resulting from all three cost functions using the 
IDVD (3rd) method are displayed in Tables 23-25, with minimum time plotted in Figures 
123-124. 
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Figure 123. Translational discrepancies for propagated IDVD (3rd) minimum time 

solution. 


1 

0.5 

0 



-0.5 


_i_i_i_i_i_i_i_ 

0 0.5 1 1.5 2 2.5 3 3.5 4 



Time (s) 


Rotational discrepancies for propagated IDVD (3rd) minimum time solution. 
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Figure 124. 


























e 

Resulting Value 

e 

Resulting Value 

e 

Resulting Value 

ei (tf) 

0.0006 

es(h) 

0.0013 

cm 

N/A 

eiitj) 

-0.0002 

^ 9 (tf) 

-0.0001 

m 

N/A 

<?3 (tf) 

0.0002 

eio {tf) 

-0.0009 

cm 

N/A 

(tf) 

-0.0009 

en (tf) 

-0.0003 

cm 

N/A 

e 5 (tf) 

0.0004 

e 12(h) 

0.0001 

m > 

N/A 

ee(tf) 

-0.0009 

ei 3 (h) 

-1.021 le-005 

cm 

N/A 

e i(tf) 

0.0055 

cm 

N/A 

H(tf.) 

N/A 


Table 23. Value of terminal conditions at the final time for IDVD (3rd) minimum time solution. 


e 

Resulting Value 

e 

Resulting Value 

e 

Resulting Value 

ei(tf) 

0.0041 

e%{tf) 

0.0003 

m 

N/A 

ei(tj) 

-0.0050 

^9 (tf) 

1.4198e-006 

cm 

N/A 

^3 (h) 

-0.0023 

e 10(h) 

-0.0020 

cm 

N/A 

£4 (tf) 

-0.0029 

en(tf) 

-0.0001 

cm 

N/A 

esitf) 

0.0021 

en(t f ) 

-0.0015 

cm 

N/A 

edtj) 

-0.0002 

ei 3 (tf) 

-0.0007 

cm 

N/A 

ei(tf) 

0.0016 

cm 

N/A 

//(h) 

N/A 


Table 24. Value of terminal conditions at the final time for IDVD (3rd) minimum energy 

solution. 


e 

Resulting Value 

e 

Resulting Value 

e 

Resulting Value 

ei(h) 

0.0019 

em 

0.0017 

cm 

N/A 

ei(tj) 

-0.0011 

em 

2.6406e-005 

cm 

N/A 

em 

-0.0002 

e 10(h) 

-0.0026 

cm 

N/A 

e 4 (tj) 

-0.0027 

en (h) 

-0.0004 

cm 

N/A 

^m 

0.0014 

e 12(h) 

-0.0005 

cm 

N/A 

em 

-0.0008 

ei3 (t f ) 

-0.0011 

cm 

N/A 

em 

0.0036 

cm 

N/A 

im 

N/A 


Table 25. Value of tenninal conditions at the final tune for IDVD (3rd) minimum fuel solution. 
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8. Summary of Close Approach Results for Baseline Maneuver 

Results for different variations of the IDVD theme pertaining to the baseline 
maneuver introduced in Chapter III and discussed throughout this dissertation are 
summarized in Table 26. 



GPOPS 

IDVD (Euler Angles) 

(vary endpoint accelerations) 







Min 

cost 






Fuel 

Final Time 

3.5086 

8.8943 

8.1182 

3.7775 

10.0000 

9.0792 

Energy 

10.4471 

0.2445 

1.1587 

5.0869 

0.3495 

0.5738 

Fuel 

20.9419 

3.6941 

2.4863 

13.3627 

4.8462 

3.9329 

CPU Time 

33 , 370.1 

86 , 040.7 

84 , 261.6 

10.6 

6.8 

11.0 



IDVD (quaternions) 

(vary endpoint accelerations) 

IDVD (quaternions) 

(vary endpoint jerk/accelerations) 


Min 

Min 

Min 

Min 


Min 

cost 

Time 

Energy 

Fuel 

Time 


Fuel 

Final Time 


8.8868 

7.5181 

3.6582 

8.9137 

8.1442 

Energy 


0.2461 

0.3875 

6.2309 

0.2458 

0.4131 

Fuel 

13.4202 

3.7177 

2.8312 

15.0982 

3.7044 

2.7552 

CPU Time (sec) 

11.5 

16.1 

32.1 

86.6 

623.0 

405.5 


Table 26. Summary of perfonnance indices and computational time for a variety of trajectory 

generation methods. 


For the cases of GPOPS and IDVD employing quaternions, Figure 125 shows the 
resulting minimum energy control histories represented in a single plot. Both the control 
histories calculated using IDVD match well with the solution provided by GPOPS. The 
similarities of the solutions are further reinforced by the small deviation in PI with the 
3rd derivative IDVD solution, IDVD (3rd), having a PI within 0.65% of the GPOPS 
solution and the IDVD (3rd-4th) solution, IDVD (3rd-4th), falling within 0.54%. 
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Figure 125. Control history of different methods for the minimum quadratic-control 

problem formulation. 


The overlaid control histories for the minimum time case are shown in Figure 126. 




Figure 126. Control history of different methods for the minimum time problem 

formulation. 
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For the minimum time case, the reduced performance of the IDVD solutions is 
mainly due to the inability to implement a bang-bang control strategy. This is because 
the IDVD results in a smooth, continuous function by construction. Although the IDVD 
(3rd-4th) strategy is more flexible, having the added capability to modify the endpoint 
jerk values, it still cannot implement the discontinuous switch of control values inherent 
to a minimum time solution for this type of problem and shows only minor improvement 
over the IDVD (3rd) strategy. Although, previous work calculating bang-bang solutions 
needed to generate an extra control of layer in order to have the required smooth control 
history for interpolation (Hurni 2009). The PI of the IDVD (3rd) solution was within 
6.5% of the GPOPS solution and the IDVD (3rd-4th) solution was within 4.3%. 

The minimum fuel solutions are shown in Figure 127. Although the IDVD 
solutions have a slightly higher PI, they avoid the rapid switching characteristics 
portrayed by the GPOPS solution inherent to the bang-off-bang control nature associated 
with the cost function described by Equations (59) and (62). The solution provided by 
IDVD is also smooth through a region that GPOPS provides a highly oscillatory solution 
incapable of being implemented. The trade off, of course, is the slightly higher PI 
associated with the solution, 13.9% increase for IDVD (3rd) and 10.8% for IDVD (3rd- 
4th). 
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Figure 127. Control history of different methods for the minimum fuel problem 

formulation. 
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VI. FURTHER SIMULATIONS AND ANALYSIS 


This chapter analyzes the rendezvous final approach problem in several different 
ways. First, the methods developed in this dissertation are applied to previously 
published examples. Concepts for closed-loop, real time operations are derived and 
demonstrated in simulation. Finally, the Inverse Dynamics in the Virtual Domain method 
is applied to a variety of scenarios with different initial conditions. 

A. COMPARISON WITH PREVIOUS STUDIES 

1. Problem Formulation in Two Dimensions 

Simulations were run to compare results with those published in previous 
literature. Specifically, the planar case of Ma (Ma et al. 2007) that is based on the 
Sakawa Shindo (SS) algorithm (Sakawa and Shindo 1980). This article attempts to solve 
the planar motion minimum time problem of a spacecraft matching position and 
orientation with a docking point that is rotating. The state vector is defined as: 

V = U y r 0 r V xr ® r f (147) 

The subscript r indicates the relative motion of the chaser with respect to the 
RSO, observed from the targets body fixed frame. There are two force controls u\ and ui, 
respectively fixed in the x and y chaser spacecraft body frame and one control torque M 3 
fixed along the z axis. The controls are bounded such that: -1 < u < 1. The initial and 
final conditions are given as: 


x r0 = (10 10 nil 1 1 1) T 
x ;/ = (1 0 0 0 0 0) T 


(148) 


Since the algorithm needs a final time, t f , as in input, the method is to solve a 

minimum fuel problem, reducing t f at each iteration. This iteration of iterations method 

can be computationally expensive, and the inability to find a solution using the SS 
algorithm does not guarantee that one does not exist. Furthennore, relative motion due to 
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the rotating reference frame (not inertially fixed) and the mean motion due to Hill’s 
equation are neglected. For more information on the SS algorithm, the reader is directed 
to (Sakawa and Shindo 1980; Sakawa 1999). Although the algorithm claims to satisfy 
Pontryagin’s Minimum Principle at every iteration, no infonnation about the costates or 
transversality conditions are supplied. The final minimum time solution arrived at by this 
method is 8.14 seconds using 200 nodes for numerical integration by the Heun method 
(Chyba, Leonard and Sontag 2000). The computation time was claimed to take several 
minutes, but not stated. Instead, the author stated that the computational time was not the 
focus of the research, only the computation of a feasible solution. 

2. Direct Method (GPOPS) Formulation and Results 

The same problem is solved by restricting the GPOPS and IDVD formulations 
derived in previous sections to two dimensions. The resulting solution for GPOPS 
provided a minimum time for the maneuver to be 7.6695 seconds. The resulting plot of 
chaser x and y coordinates fixed are presented in the RSO body frame, as in Ma (2007), is 
shown in Figure 128. The control histories are shown in Figure 129, as well as the 
switching conditions developed in Chapter III.B. The controls obey the bang-bang nature 
stated in Ma (2007) and Chapter III of this dissertation, which are dictated by the 
switching functions. The controls presented in Ma are said to approach the bang-bang 
structure but appear rather smooth. Figures 130 and 131 show the time history of the 
endpoint conditions, the difference in position and velocity of the docking points, as well 
as angular rate and orientation of the vehicles, converging to zero. The transversality 
conditions and Hamiltonian shown in Figures 132 and 133 further reinforce the optimal 
nature of the control for the minimum time cost. Still, the major drawback of this 
formulation and method is that the calculation time of the solution was 1,416.2 seconds 
(23.6 minutes) for only a 50 node solution (not shown), and the calculated solution for 
200 nodes, shown in Figures 128-133, required 40,659.0 seconds of computational time. 
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4 6 8 10 12 14 

X coord RSO Body Frame 


Figure 128. 2D Case (GPOPS): optimal minimum time trajectory of the chaser in the RSO 

body frame. 



Time (s) 


Figure 129. 2D Case (GPOPS): optimal minimum time control history. 
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Figure 130. 2D Case (GPOPS): history of the translational endpoint conditions. 
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Figure 131. 2D Case (GPOPS): history of the attitude endpoint conditions. 
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Figure 132. 2D Case (GPOPS): history of the transvervsality conditions. 
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Figure 133. 2D Case (GPOPS): history of the Hamiltonian. 


For comparison purposes, the GPOPS method was run with a minimum fuel cost 
based on Equation (93) but with a multiplier of 0.005 (Ma et al. 2007) and a maximum 
final time of 8.14 seconds. The solution took 1,953.1 seconds (32.6 minutes) to calculate 
a 50 node solution and provided a PI of 0.0574 (compared to the 0.0788 PI result from 
the SS based method). 
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3. IDVD Formulation and Result 

The IDVD method was also used to calculate the trajectory based, shown in 
Figure 134, on the 2D problem posed by Ma (Ma et al. 2007). The associated control 
history is shown in Figure 135. The endpoint conditions are shown in Figures 136 and 
137. Again, the same formulation and methodology was used as in Chapter V.A, varying 
acceleration at the endpoints, but the problem was constrained to planar motion. The 
resulting costs and plots were generated using 200 points (the resulting IDVD is analytic, 
so any number of nodes or waypoints can be used to present the final solution). This 
final time calculated, 8.7525 seconds, is higher than both the previous methods, but the 
solution only took 10.3 seconds to compute. When the IDVD (3rd-4th) method is used, 
the final time to complete the maneuver is decreased to 8.2086 seconds, but the 
computational time to arrive at the solution increase to 214.7 seconds. 
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Figure 134. 2D Case (IDVD): minimum time trajectory of the chaser in the RSO body 

frame. 
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Figure 135. 2D Case (IDVD 3rd): minimum time control history. 
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Figure 136. 2D Case (IDVD 3rd): history of the translational endpoint conditions. 
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Figure 137. 2D Case (IDVD 3rd): history of the attitude endpoint conditions. 

The increased flexibility of the IDVD (3rd-4th) basis functions permits the 
solution to behave more like the bang-bang structure of the optimal solution, allowing for 
a reduced final time for the maneuver to complete. 

4. Comparisons of Trajectory Generation Methods for 2D Example 

Table 27 shows a comparison between methods for the minimum time solution. 
200 nodes were chosen for the calculation and a factor of 0.005 was multiplied by the 
minimum fuel cost to be consistent with previous literature (Ma 2007). Although IDVD 
(3rd) has a greater final maneuver time, the computational time is drastically reduced. 
The computational time, simply stated as several minutes, and final minimum energy cost 
of the SS solution presented by Ma (2007), were not reported. The solution providing the 
best minimum time cost was the bang-bang solution solved for by GPOPS, but the 
solution came at an extensive computational price. 
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GPOPS 

IDVD 

(3rd) 

IDVD 
(3rd-4th) 

Ma (SS) 

cost 

Min 

Time 

Min 

Time 

Min 

Time 

Min 

Time 

Final Time 

7.6695 

8.7525 

8.2086 

8.1400 

Energy 

11.3266 

5.9307 

6.9566 

N/A 

Fuel* 

0.1135 

0.0779 

0.0824 

0.0788 

CPU Time (sec) 

40,659.0 

10.3 

214.7 

N/A 


Table 27. Summary of performance indices and computational time for a 2D scenario. *A factor 
of 0.005 was multiplied by the minimum fuel cost to be consistent with the fonnulation 

presented by Ma (2007). 

B. CYCLICAL NATURE OF PROBLEM SOLUTION PERFORMANCE 
INDEX 

The performance index of rendezvous of a spacecraft with a tumbling object 
exhibits a cyclical nature as seen by Figures 138 and 139 for two scenarios. Figure 138 
shows the cyclic nature of the PI with respect to a time value associated with waiting to 
start the maneuver. If the initial angular velocity of the RSO occurs around a principal 
moment of inertia, then the motion of the docking point is periodic, containing only 
circular motion in the plane perpendicular to the angular velocity vector. A simple 
example of fixing X =1 for the entire maneuver and then varying the wait time before 
maneuver start (in essence is the same as varying the initial orientation of the RSO), 
shows that the PI associated with completing the maneuver is periodic with time. Figure 
139 shows the cyclic cost nature for the case of the docking point motion, not being 
constrained to a plane perpendicular to the initial angular velocity vector of the RSO. 
Such is the case of the RSO having a nonidentity inertia matrix and initial angular 
velocity vector that is not coincident with a principal moment of inertia, the PI is cyclic, 
but does not repeat. For this case, the time value shown is the total time of the maneuver. 
For both cases, there is an inherent cyclical nature to the PI pertaining either the wait time 
until commencing the maneuver or the total time of the maneuver. This added 
complexity demonstrates the increased influence of the maximum final time on the 
solution. For any given circumstance, allowing more time to perform the maneuver may 
not necessarily decrease the PI for minimum fuel or minimum energy maneuvers as 
would be the case if rendezvous was to a fixed point in space. For this reason, initial 
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guesses should be chosen judiciously and when possible, as stated in Chapter V.B, and 
should provide a feasible solution. This way, since the kinematic and dynamic equations 
are always satisfied though IDVD, there is always a solution to implement at any given 
time should the user want to terminate the optimization routine. 



Figure 138. Energy and fuel costs for rendezvous with a tumbling RSO with symmetric 
inertia matrix. The time shown is a wait time until the start of the maneuver. 



Figure 139. Energy and fuel costs for rendezvous with a tumbling RSO with an asymmetric 

inertia matrix. The time shown is the total maneuver time. 
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C. FLEXIBILITY AND ROBUSTNESS OF IDVD TECHNIQUE 

The purpose of this section is not to generate an exhaustive list of potential 
advantages and implementations of the IDVD methods described in the previous section, 
but to further educate the reader on concepts in order to facilitate new ideas. The 
minimum energy example was selected for further analysis and implementation based on 
its nature to minimize control effort over the entire maneuver. It has inherent fuel¬ 
minimizing qualities as well as other attractive attributes. The minimum energy cost 
based maneuver does not have the trait of commanding maximum (or saturated) controls 
during the tenninal phase of docking. The extra control margin makes it the most 
responsive and safe if sudden changes in trajectory are needed, a potential necessity if 
closed loop implementation is to be considered. Having low control action also lowers 
potential effects of plume impingement as opposed to maximum thrust. These combined 
attributes make it the preferred choice for implementation. 

1. Real Time Trajectory Reshaping 

Because the IDVD solution is analytic in nature and incorporates prescription of 
the endpoint values as well as their derivatives, it allows for the possibility of 
implementation in a closed-loop fashion. If, for example, the minimum energy trajectory 
(which does not saturate the actuators at any point) solution was calculated using the 
IDVD method, the resulting varied parameters would vary the higher level derivatives at 
the endpoints, since the state values and their respective first derivatives were set based 
on the projected state of the RSO at the final time as well as the initial conditions of the 
chaser craft. Even after a solution is calculated, if updated information about the current 
states of the RSO were obtained and the final state values of the RSO were to change 
(due to unmodeled dynamics or disturbances), this infonnation can be reintroduced into 
the trajectory generation method to ensure that the endpoints of the trajectory met the 
desired conditions. Mathematically speaking, the varied parameters, already solved for, 
would not change, but, keeping the conditions at the beginning of the trajectory fixed, the 
specified conditions at the ending point, tf, would be tweaked in order to match the best 
projected conditions of the RSO docking point, resulting in new coefficients for the 
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polynomials and therefore instantly reshaping the trajectory. The new trajectory may not 
be optimized with respect to the slight changes in state variables, but it would be the best 
solution based on the given constraints (having a feasible trajectory that ends with the 
desired conditions on the chaser and its docking point) without going through another 
iteration process. This is in contrast to current iterative real time closed loop optimal- 
control studies that are performed in simulation, (McFarland et al. 2009), and assume 
solutions can be obtained in whatever specified update rate is needed, stopping the 
simulation to compute the optimal control, which can take several minutes to hours 
depending on computer perfonnance and desired solution resolution. The following two 
examples are used to highlight the noniterative reshaping trait of the proposed IDVD 
controller. In the first scenario, the inertia matrix of the RSO is thought to be known as I 
= diag([3,1,2]); therefore, it is used to project the states of the RSO at tf based on the 
initial conditions of the RSO and Euler’s equations for rotational motion. For this 
exercise, we assume the solution of the varied parameters that optimize the rendezvous 
maneuver is already known and the conditions are such that the problem is the same as 
discussed in previous sections based on Table 1. The only difference is that now the 
actual inertia of the RSO is I = diag ([1,2,3]). The second scenario involves correct 
knowledge of the RSO inertia, but the presence of an unknown constant thrust in the 
body frame (such as a thruster that is failed on). New (ideal) RSO state infonnation is 
available at a rate of 10 Hz, allowing the trajectory generator to reshape in real time 
based on the new projected endpoint values. A block diagram of the process is shown in 
Figure 140. 

In the previous sections, the IDVD Solver coupled with the Trajectory Generator 
was treated as the same system. It should be noted that, in fact, they can be separated in 
order to exploit specific advantages to the IDVD method. The key difference between 
the IDVD Solver and the Trajectory Generator is the IDVD Solver iterates on the varied 
parameters (the higher level derivatives of a given trajectory at the endpoints) and uses 
the Trajectory Generator to solve for the polynomial coefficients and generate the spatial 
and time trajectories of the system states (as well as the controls). From there, a PI can 
be associated with the recently calculated trajectory and the IDVD Solver can iterate on 
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those varied parameters to reduce that PI. In summary, the entire IDVD makes calls to 
the Trajectory Generator. The Trajectory Generator, when examined alone, takes the 
endpoint conditions (position and velocity) and higher order derivatives (eg. acceleration) 
for attitude and translation and solves for the polynomial coefficients using Equations 
(121)—(128) for attitude and Equations (136)—(139) for translation. The idea being that 
once the IDVD Solver has values for the varied parameters, the most current information 
on the RSO docking point can be used to reshape the trajectory, ensuring endpoint 
conditions are met. Furthermore, it can be implemented in such a manner that also 
provides the controls necessary to track this trajectory, using the inner-loop controller to 
take care of small errors (Yakimenko et al. 2008). 



Figure 140. Block diagram of real time trajectory reshaping concepts using IDVD. 

a. Solution Results with Inertia Uncertainty 

The trajectory resulting from the closed-loop implementation is shown in 
Figure 142. For this case, the trajectory generator has incorrect information about the 
inertia of the RSO, specifically it is using a diag ([3,1,2]) inertia matrix (the same as 
previous examples) when the actual inertia matrix is diag ([1,2,3]). This is a more 
extreme case of misidentification and was chosen to best illustrate the concept. In 
actuality, the inertia matrix information of the RSO employed by the trajectory generator 
would most likely be more accurate, or simply use an identity inertia matrix. 
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In Figure 141, the actual trajectory of the RSO docking point is illustrated 
by a dotted red line with a circle marker highlighting a subset of points at which updates 
are provided to the RSO. The dash/dot lines illustrate some of the projected trajectories 
of the RSO docking point based on the most current angular position and velocity 
information, propagated with the (incorrect) RSO inertia matrix currently used by the 
trajectory generator. The dotted line shows the projected endpoint of the RSO docking 
point as it progresses though each iteration. The projected RSO docking point and the 
actual RSO docking point converge to the same value as at tf, since the amount of time 
for the trajectory generator to project into the future, based on incorrect inertia 
information, where the RSO docking point will be decreases. The thin colored lines 
show a subset of the reshaped trajectories at each timestep. The blue dots show the 
resulting overall trajectory for the chaser vehicle, overlaid with intermittent models 
showing the resulting attitude of the vehicle. 
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Figure 141. Results of the final trajectory for the close approach example having a misidentified inertia matrix IDVD real time reshaping. 
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Figure 142 shows a close-up, expanded view of the previously described trajectories. 
Information shown is the actual trajectory (red, circle), the projected endpoint (black, 
dash) and a sample of the overall projected trajectories based on state updates (green, 
dash/dot). Figure 143 illustrates the evolution of the current trajectory over time and the 
position of the chaser CM. 
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Figure 142. Exploded view of the RSO docking point information for a misidentilicd 

inertia matrix example. 
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Figure 143. Spatial views of the resulting trajectories for the close approach example 
having a misidentified inertia matrix IDVD real time reshaping. 

b. Solution Results with Unknown Constant Torque 

For the next example, it is assumed the inertia matrix of the RSO is known 
by the trajectory generator, but an u nkn own torque exists, acting on the vehicle in the x- 
body frame. This can be thought of as a thruster that is stuck in the On mode, thus 
generating the disturbance. The resulting plots, similar to those in the previous section 
are shown in Figures 144 and 145. Figure 144 shows several of the reshaped trajectories, 
along with the overall trajectory for the chaser vehicle, while Figure 145 shows the time 
history for each translational component of the Chaser vehicle, as well as the evolution of 
the reshaped trajectories. As in the previous example, the projected and actual endpoint 
of the RSO states converge to zero as the chaser vehicles converges to the final state. 
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Figure 144. View of the close approach trajectory example with an unidentified, constant torque using IDVD real time reshaping. 
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Figure 145. Spatial views of the resulting trajectories for the close approach example 
having an unidentified constant torque using IDVD real time reshaping. 

2. Resolving and Reoptimizing the Problem in Real Time 

If the solution of the IDVD could be computed in real time, with the definition of 
real time based on the requirements of the user/customer and the estimation attributes of 
the sensor system, the problem can be recomputed as the trajectory is traversed. An 
illustrative example of this concept is given in Figure 146. In this case, a switch is 
installed that allows the initial conditions to be updated based on the most current system 
states. The problem can then be resolved using IDVD and the new resulting trajectory 
immediately updated. The key difference between this implementation and the reshaping 
approach is that the trajectory is resolved and reoptimized based on the most current 
states of the RSO and the chaser vehicle. Even though for this case the trajectory is 
resolved, one of the greatest benefits between the IDVD method is still that if the 
endpoint conditions change, the IDVD can use this information to tweak the trajectory. 
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Other methods, specifically pseudospectral, must resolve the entire problem in order to 
obtain a trajectory that, when the controls are propagated, would finish the maneuver in 
the correct position with the desired conditions. 
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Figure 146. Block diagram of real time trajectory recalculation concepts using IDVD. 

This concept is demonstrated on the same example stated in Chapter VII.D.l.b, 
employing the inertia uncertainty. The trajectory generator is hampered by not knowing 
the true inertia of the RSO, but is receiving state updates and recalculating the trajectory 
at a rate of 1Hz. This allows the current trajectory to be based on the most current kn own 
states of the chaser and the RSO. The switch allows for the trajectory generation method 
to use the rapid-reshaping concept described in Chapter VII.D.l with a 10 Hz update rate 
for the final 1 second, or endgame phase, of the trajectory since uncertainties in guidance 
may need to be corrected much more rapidly (Vaidyanathan et al. 2001). This also 
demonstrates the IDVD’s ability to be integrated with other guidance systems to 
construct the best overall system. 

Figure 147 shows the plots of each recalculation (top), as well as the overall 
resulting trajectory (Bottom). The multicolored segmented lines in the top figure 
illustrate each trajectory recalculation, occurring at 1-second intervals shown by the blue 
dots. The segmented line on the lm sphere shows black “x” marks representing the 
current projected location of the docking point at the calculated tf for each optimized 
trajectory. 
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Figure 147. (Top) Overlaid recomputed trajectories. (Bottom) Resulting overall trajectory 

based on the most current recalculation. 


3. Specific Waypoint Spacing in the Time Domain 

Another advantage to the specific IDVD formulation, and an enabling attribute 
that allowed the implementation of the closed-loop implementation concepts from the 


175 










































































previous sections, is the freedom to choose node or waypoints. In previous work 
(Bevilacqua, Romano and Yakimenko 2009), integrating the inverse of speed factor that 
was based on switching times led to the user to be at the mercy of node spacing that is 
dependent upon the history of the speed factor itself. This could lead to gaps in control 
and trajectory information that make the solution incapable of implementation. This is 
also seen in application of pseudospectral direct methods (Humi 2009) while using a 
commercial program for direct optimization of trajectories and having to add layer of 
control (and complexity) in order to correctly interpolate control information in between 
the node points of the solution. Because of this, for the problem studied concerning a 
wheeled rover-like vehicle, the user is forced to use velocity and heading as a control 
vector, a layer above the intuitive controls of acceleration and angular rate. An example 
where the constraints on the previous problem are relaxed to illustrate how the nodes of 
the solution get bunched together in regions where the speed factor is the greatest is 
considered. This is because the node spacing is uniform with respect to the virtual 
argument t , as shown in the bottom plot of Figure 148. 
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Figure 148. Control history of a representative trajectory with uniform node spacing in the 

virtual domain. 
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The question posed is what dispersion of node points in the virtual domain will 
guarantee uniform node spacing in the time domain. For the methods described in 
Chapters III and IV, the speed factor was chosen such it has an analytic expression for the 
integral of its inverse, which maps a specific value of t to the time domain as shown 
below: 


t = 



(149) 


This results in the ability to express z as a function of t and vice versa. 

Once the final time is obtained from the initial solution, a set of evenly space 
nodes (or arbitrarily spaced nodes based on any inputs from the user) set in the time 
domain can be instantaneously converted to the virtual domain. This new set of nodes in 
the virtual domain, when implemented, leads to the evenly spaced nodes in the time 
domain as shown in Figure 149. 
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Figure 149. Control history of a representative trajectory with uniform node spacing in the 

time domain. 
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D. EXTENDED EXAMPLES USING INVERSE DYNAMICS IN THE 
VIRTUAL DOMAIN 

Up to this point, the same initial conditions have been used to analyze the optimal 
solution to a close approach problem, compare with a suboptimal rapid-trajectory 
generation technique, and demonstrate potential real time close-loop applications. The 
IDVD (3rd) and IDVD (3rd-4th) method generated in the previous chapter is applied to 
several pseudorandom initial chaser conditions and scenarios. First, 10 samples were 
developed to apply the IDVD technique. From these 10 samples, the initial velocities of 
the spacecraft were assumed to be zero and the initial orientation of the chaser is assumed 
to be coincident with the orbit frame; both assumptions were taken from initial conditions 
used by McCamish (2007). The pseudorandom initial positions of the chaser are chosen 
by first detennining a vector them multiplying it by a uniformly distributed random 
number between 5 m and 10 m as shown in Equation (150). 

x = - rand x 2 )cos (rand e )rand r 

v = ^j( 1 - rand x 2 j sin (rand g )rand r (150) 

z = rand/and r 

with rand x e [0 1] , randg e.\ 0 2n), rand r e [5 10] having unifonn random 

distributions. The desire was to have an equally distributed set of points in range (from 
5-10 m) and all directions, understanding that based on the given formula, the points 
would not necessarily be uniformly distributed about the volume housing those points. 
The solutions to the first ten generated initial conditions using IDVD (3rd) are illustrated 
in Figure 150. Again, the initial conditions of the RSO were kept consistent with the 
previous examples from this manuscript, having a I = diag[3 1 2] and initial angular 
velocity of 0.25 radians/second in the y and z body frame and initial q = [0 0 0 1]. For 
the first set of ten shown in Figure 150, the average computational time was 21.0 
seconds. 
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Figure 150. Trajectories for sample rendezvous problem for the first 10 random starting 

chaser positions with zero initial relative velocity. 

Next, over 1000 pseudorandom initial conditions were generated with varying 
initial relative velocity as well as position. The maximum magnitude of the initial 
velocity was limited to 1 m/s to be consistent with “relative speed limit” constraints 
placed during close proximity operation again in the complementary work of McCamish 
(2007) as well as having the body frame aligned with the orbital frame. This would be 
consistent with a GNC concept utilizing McCamish’s Artificial Potential Function 
modified algorithm up until it is time to perform the close approach for docking 
maneuver when the RSO has nonzero angular rates. The solutions to the first 10 samples 
are shown in Figure 151. Notice the presence of initial velocity for the chaser vehicles 
results in a modified path to the final conditions. This is due to the fact the resulting 
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polynomial trajectory must shape itself to have a directional component at the beginning 
of the maneuver that will accommodate the different initial conditions. 


X" f 



Figure 151. Trajectories for sample rendezvous problem for the first 10 random starting 

chaser positions with nonzero initial relative velocity. 

The complete distribution of initial conditions is shown in Figures 152-155. The 
maximum allowable final time was also increased to 15 seconds for the ensuing 
calculations. 
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Figure 152. Distribution of the 1000 pseudorandom initial conditions tested. 
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Figure 153. Planar view of the 1000 pseudorandom initial conditions tested. 
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Figure 154. 


Figure 155. 
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Magnitude of Initial Range from RSO (m) 


Distribution of the initial chaser range for the 1000 pseudorandom samples 

generated. 
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Magnitude of Initial Velocity (m/s) 


Distribution of the initial chaser velocity for the 1000 pseudorandom samples 

generated. 
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The IDVD (3rd) method was run using the 1000 sample and the mean time to calculate 
the solution was 34.3 seconds. The mean time to complete the maneuver was 14.3081 
seconds and the mean energy PI was 0.2289. The same 1000 initial condition set was 
also run using the IDVD (3rd-4th) method. The mean of the energy PI was 0.2112, but 
the average run time was 228.9097 seconds. A summary of the results is shown in Table 
28 for initial conditions on relative position and velocity pertaining to Figures 154 and 
155 and Table 29 pertaining to initial relative positions based on Figure 154 but starting 
with zero initial relative velocity. Table 30 shows the results for varying the initial 
conditions on position and velocity of the chaser vehicle according to Figures 154 and 
155 but also varying the initial angular velocity of the RSO based on Figure 156. 
Although the IDVD (3rd-4th) method results in an overall average cost over the sample 
of initial conditions tested, the computational time was substantially more. 



IDVD 

(3rd) 

IDVD 

(3rd-4th) 

mean values 

Min 

Energy 

Min 

Energy 

Final Time 

14.3081 

14.5191 

Energy 

0.2289 

0.2112 

Fuel 

4.3124 

4.2516 

CPU Time (sec) 

34.3 

228.9 


Table 28. Summary of perfonnance indices and computational time for IDVD (3rd) and IDVD 
(3rd-4th) method using 1000 pseudorandom samples for varying the initial conditions 

of the chaser position and velocity. 



IDVD 

(3rd) 

IDVD 

(3rd-4th) 

mean values 

Min 

Energy 

Min 

Energy 

Final Time 

14.3094 

14.6107 

Energy 

0.1742 

0.1631 

Fuel 

3.8773 

3.8544 

CPU Time (sec) 

32.9 

191.0 


Table 29. Summary of perfonnance indices and computational time for IDVD (3rd) and IDVD 
(3rd-4th) method using 1000 pseudorandom samples for varying the initial conditions 
of the chaser position but having the initial relative velocity be equal to zero. 
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Magnitude of Initial RSO Angular Velocity (rad/s) 

Figure 156. Distribution of the initial RSO angular velocity for the 1000 pseudorandom 

samples generated. 



IDVD 

(3rd) 

IDVD 

(3rd-4th) 

mean values 

Min 

Energy 

Min 

Energy 

Final Time 

14.9716 

14.9799 

Energy 

0.2138 

0.1689854 

Fuel 

3.3867 

3.2472 

CPU Time (sec) 

20.0 

228.5 


Table 30. Summary of performance indices and computational time for IDVD (3rd) and IDVD 

(3rd-4th) method using 1000 pseudorandom samples for varying the initial conditions 
for position and velocity of the chaser and initial angular velocity of the RSO. 
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VII. CONCLUSION 


This section summarizes the contributions and findings of this dissertation and is 
organized as follows: First, a summary of the research is presented. This is followed by a 
summary of the contributions of this research followed by a discussion of the advantages 
and disadvantages pertaining to the different methods. Finally, a recommendation for 
future study is presented. 

A. SUMMARY OF RESEARCH 

A six DoF, 20-state model of two spacecraft rendezvous is developed, one of 
which was controlled, the other considered to be passively tumbling. A solution is 
obtained for the problem of close approach, up to the point of contact, using a direct 
optimal control method. The solution is then verified as optimal by way of an indirect 
method based on the MP. Next, a trajectory generation method for spacecraft 
reorientation is developed, based on a quaternion construction of IDVD. This new 
construction enables implementation of an IDVD trajectory generation method for the 
problem of a spacecraft performing a close approach maneuver to a tumbling object. 
Finally, the advantages of the IDVD method were exploited and demonstrated in 
simulated scenarios that employ closed-loop feedback. 

B. SUMMARY OF CONTRIBUTIONS 

Contributions of this dissertation are considered to have two distinct categories of 
impact, theoretical contributions, enhancing the overall understanding contribute to the 
overall body of knowledge, and practical contributions, advancing the current state of the 
art and moving the technology closer to implementation. They are both summarized in 
the next several paragraphs. 

Summarizing, the minimum quadratic-control, minimum fuel and minimum time 
continuous optimal control 3D spacecraft rendezvous-docking problems are fonnulated 
for the first time in literature, and addressed using a direct collocation method, GPOPS. 
For both problems, the desired optimal trajectory of chaser spacecraft with respect to a 
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tumbling RSO is sought such that the desired docking points match in position and 
velocity. Moreover, the solutions obtained were verified in the numerical simulations 
based on the Minimum Principle. Those solutions included derivation of the adjoint 
equations, formulation of the necessary conditions for the optimal solution, and synthesis 
of the optimal controls. The results obtained from using the direct collocation method— 
in this case the open source software of GPOPS—are very close to those obtained using 
with the Minimum Principle. It was also found that path constraints are necessary when 
solving for the optimal trajectory in order to prevent undesired collision of spacecraft. 
This was shown by the active path constraint that results in discontinuous costates upon 
contact with the constraint boundary. 

As expected, the GPOPS direct collocation (pseudospectral) method was found to 
be quite reliable and yielded the results relatively fast. Moreover, these results were also 
validated via propagation of the states based on the system dynamics, and optimal 
controls were found. That is what would be needed for possible onboard implementation 
of the developed algorithms. However, the GPOPS direct collocation method proved to 
be unable to produce the results fast enough so that it could be used in real time. Hence, 
in the nearest future, it only leaves off-line open-loop option to be possibly used on the 
real spacecraft. But even then, it was found that in case of a singular control, the method 
returns somewhat infeasible results that cannot be implemented. The main simplifying 
hypothesis, considered in the numerical simulation presented in this section of the 
research, was the spherical inertial symmetry of both chaser and target spacecraft, as well 
as thrust control limits being applied in the orbital frame regardless of chaser spacecraft 
orientation. In order to achieve feasible controls in real time, this assumption needs to be 
removed, along with investigating the usage of different combined Pis and other direct 
methods based on inversing of the dynamics of the problem. 

This research then continued work on spacecraft rendezvous by modifying the 
previous problem formulation to include body mounted translational thrusters on the 
chaser spacecraft, as well as extensions on the final conditions. These new final 
conditions involved matching the velocity and position of the spacecraft docking points, 
as well as matching angular velocity and orientation. Inertia parameters on the RSO were 
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also set to different values, removing the assumption of an identity inertia matrix, in order 
to stimulate dynamic angular rates in all three axes. The problem was addressed again 
using the most plausible pseudospectral methods. In addition, the optimal control 
structure was synthesized and analytical expressions for the transversality conditions 
were derived (exploiting the Minimum Principle) to compare those with the calculated 
solution. In this section of research, three Pis were examined. It was found that the 
minimum time solution exhibits the expected bang-bang control and the minimum fuel 
solution exhibits bang-off-bang characteristics. Assuming controls in the body frame as 
opposed to the orbital frame (as in the previous work by the authors), while making the 
problem more complex, seems to alleviate the singular control that appeared in the z- 
translation, and therefore, made the GPOPS solution feasible. Another finding was that 
instead of addressing the minimum fuel problem, it was wiser to consider a minimum- 
energy problem instead. The reason for this is that although the minimum fuel solution 
happened to be about 30% more effective, it would be difficult to implement the control 
in practice because of the bang-off-bang control structure. The minimum energy solution, 
which lasts longer (but follows a similar path), assures smooth, easy-to-follow control 
time histories. Furthermore, the minimum energy solution does not involve having the 
actuators saturated at the final time (a characteristic of bang-off-bang maneuvers for 
minimum fuel or bang-bang for minimum time), increasing the safety of the maneuver by 
having control margin available if a rapid corrective maneuver is needed and reducing the 
likelihood for or large disturbances due to plume impingement. The minimum time 
solution (following a different path) does not contact any state constraints and therefore 
has continuous costates, while the minimum-control and minimum energy solutions 
contact the keep-out zone of the RSO, while approaching for docking. Finally, it was 
shown that the required CPU time greatly exceeds that of the maneuver itself, which as of 
now excludes pseudospectral direct collocation methods from being a candidate for 
onboard application. Instead, exploitation of other approaches, willing to sacrifice a 
fraction of the PI but decrease the required CPU time by several orders of magnitude, 
were investigated, as summarized below. 


187 



In the case of reorientation maneuvers, the IDVD method with the novel 
quaternion approximation functions proposed in this dissertation allows computing 
feasible solutions fast enough to be used onboard satellites for on-line computation of 
slew maneuvers. Moreover, because of the smooth controls' histories it can be 
implemented in the control schemes involving a feed-forward loop. For this specific 
implementation a form of the speed factor was chosen based on a combination of 
nonlinear equations as opposed to integrating dynamic equations based on having 
switching times of the controls act as varied parameters (Yakimenko 2000). Because of 
this, compared to the true time-optimal solutions, the IDVD trajectories do not have the 
capability to generate a bang-bang control solution, which results in a slightly worse PI. 
However, smooth controls benefit other mission preferences of having desired rates at the 
endpoints as well as guaranteeing the endpoint constraints of the maneuver are always 
satisfied. In addition, it is a definite advantage in rapidly changing acquisition or tracking 
scenarios and when the slewing spacecraft possesses low frequency flexible modes. This 
formulation is later applied to situations where the attitude is coupled with other 
dynamics such as translational motion in rendezvous and docking applications. In this 
case, a simple eigenaxis slew would not meet mission criteria such as matching rotational 
motion. 

IDVD is further applied to the maneuver that has both translational and attitude 
motion, such as the case with the focus of this research involving close approach to a 
tumbling object. The IDVD method developed is highly flexible, allowing for matching 
any mixed set of endpoint conditions, generating user specified waypoint spacing, and 
providing a mechanism for real time closed-loop trajectory reshaping. Based on the 
trajectory generation results for the benchmark 3D scenario, the minimum-control 
solution was chosen for further study and implementation. This was due to its inherent 
characteristic to conserve fuel, while providing specific safety qualities, such as not 
commanding actuator saturation (which would limit its ability to do rapid reactive 
maneuvers in case of closed loop implementation of higher-level emergency abort). 
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Within the IDVD method, there are still variations in trajectory parameterization 
techniques. The IDVD 3rd method employs varying the 3rd derivative at the endpoints 
of the trajectory while the IDVD 3rd-4th method uses both the 3rd and 4th derivatives as 
varied parameters. Because of this, the IDVD 3rd-4th method can provide a final 
trajectory that is lower in PI but at the cost of computational time. The IDVD 3rd method 
can obtain a solution that is within 10% of the IDVD 3rd-4th method while only using 
15% of the computational time, based on results from 1000 samples with random initial 
conditions. 

C. CONSIDERATIONS AND RECOMMENDATIONS 

Trajectory generation by direct collocation methods have the benefit of providing 
a mathematically verifiable optimal solution. The major drawback is the excessive 
computational time needed on current computing systems. The resolution of the solution, 
and computational time, is highly dependent on the amount of nodes used to solve the 
problem. Furthermore, the user cannot specify where the nodes are placed and this leads 
to complex interpolation schemes or adding an extra layer of control, increasing the 
complexity of the problem formulation in order to have sufficient control infonnation to 
implement. The benefits seen by IDVD methods are the rapid computational time that 
allows for a feasible solution to be generated, potentially onboard a spacecraft and in 
closed-loop. Because of the analytic nature of its solution, high node trajectories can be 
generated without sacrificing computational time. The ability to reshape itself when 
provided updated information on the RSO states and the docking point reinforce the 
overall robustness of the IDVD method for trajectory planning. By increasing the order 
of the basis polynomials, thereby, increasing the varied parameters, solutions can be 
obtained that better minimize the PI, but at the cost of higher computational times. 
Furthermore, the node spacing of the control solution can be completely specified by the 
user, therefore, eliminating the previously discussed error incurred by interpolation of 
potentially complex control profiles or adding extra layers of control because of the lack 
of flexibility in the node spacing in the control solution set by the direct collocation 
methods, which are not guaranteed to be continuous. It should be noted that using the 

IDVD formulation, constraints can be set on higher derivatives than the level of the 
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control vector. For example, constraints on the jerk profiles can be enforced, while still 
using the acceleration profile as the control vector. This attribute, which is in complete 
contrast other direct methods, is based on the analytical formulation of the trajectory. For 
example, to set constraints on jerk using direct collocation method techniques, the control 
vector would have to at least be based on jerk (or even higher derivative) and not 
acceleration, again adding additional levels of complexity. 

Drawbacks to the IDVD method are that, while they are optimized to give the best 
trajectory based on the constraints of the basis polynomials, they cannot give the truly 
optimal solution described by the infinite dimensional optimal control problem, nor 
would an optimal solution be verifiable through Pontryagin’s MP due to the lack of 
costate information. Also, maneuver duration is an issue that currently plagues both the 
direct and IDVD methods. For IDVD this is due to the additional coupling of the attitude 
and translational motion through the speed factor. If the speed factor is made small to 
accommodate constraints on translational motion and there is a nonzero requirement on a 
higher order derivative at an attitude maneuver endpoint (note this issue does not exist for 
translational maneuvers or rest-to-rest maneuvers), then the attitude maneuver may 
"wrap" or rotate beyond 2n in order to perfonn the maneuver to meet the endpoint 
requirements. This is because, for the given set higher order derivative in the time 
domain, a reduction of speed factor at the end point (which may be needed to 
accommodate maneuvers with excessive times) would lead to a larger value for the 
respective higher order derivative in the virtual time domain since X appears in the 
denominator of the conversion multiplier from time to virtual domains. This also 
translates into a performing a larger maneuver in quaternion space to accommodate the 
higher derivative value in the virtual domain. While the maneuver would still be 
completely feasible, the extra rotation may be undesired in the final approach. Even 
though this is not an issue for the implementation discussed in this dissertation, it would 
need to be addressed for application on a subset of missions that have drastically different 
timescales associated with the attitude and translational maneuver. Although this is 
easily rectified by simply setting constraints on the speed factor so the extra rotation does 
not happen, it may hinder the resulting PI. Another approach is employing a trajectory 
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generation scheme for maneuvers starting farther away (far approach), which simply 
holds attitude fixed in the orbital frame (McCamish 2007) or pointing to the RSO, then 
switching to the final approach scheme discussed in this dissertation. This is in line with 
operational methodology that would want to be pointing at the target as long as possible 
in order to get the most up-to-date information on the tumbling RSO before performing 
the rapid final approach maneuver at the last possible instant. On the other hand, for the 
direct collocation methods, maneuvers that are particularly long in length are subject to 
large time gaps in between nodes, where no state or control infonnation is available. 
Extra caution needs to be taken when employing such a solution. 

Summarizing, two novel guidance strategies for final approach to a tumbling 
object are developed in this dissertation. One is based on Inverse Dynamics in the 
Virtual Domain, employing a newly fonnulated quaternion approach, and the other is a 
novel formulation that employs direct collocation methods for optimal trajectory 
generation. While the direct collocation method formulation is too computationally 
expensive to be employed, it provides a mathematically verifiable optimal baseline for 
the maneuvers. While the method based on IDVD cannot exactly match the optimal 
solution, it can provide a feasible that is optimized stemming from set of basis functions 
with a fraction of the computational time of the direct method. This and other attractive 
features discussed, allows for exploitation of several real time implementation concepts. 
It is recommended that the IDVD method be employed for situations where the 
computational cost of the solution outweighs the PI associated with the state and control 
histories. The direct collocation, or pseudospectral, method should be used in 
circumstances where computational time is not a concern, but a truly optimal solution 
based on the state and control variables is desired, such as base lining scenarios and 
finding limits of specific technologies. 
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D. POSSIBLE FUTURE DEVELOPMENTS 

The following research problems remain open for possible follow-up efforts: 

1. Investigate effects of different sets of basis function for trajectory 
generation and overall effect of using IDVD methods to seed initial 
guesses for pseudospectral optimal control solvers. 

2. Analyze performance of trajectory tracking in simulation. 

3. Experimental validation on real autonomous vehicles. 

4. Integration with current research involving navigation and target 
identification. 
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APPENDIX. QUATERNION PROPERTIES 


A quaternion is a four-dimensional vector that is used to express orientation of a 
rigid body. It is composed of a scalar, 5 , and a vector portion, v. There are several 
representations and formulation of quaternions and their associated properties in 
literature. Throughout this dissertation, the properties in this Appendix will be employed 
(Titterton and Weston 1997; Sarkka 2007). 


Quaternion Expression: 





UJ 


(151) 


Quaternion Length: 




(152) 


Quaternion Logarithm: 



(153) 


Quaternion Exponential: 


exp(q) = exp(s) 



(154) 
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Quaternion Conjugate: 




/ 

-q 2 

-q 3 

) 


Quaternion Inverse: 



Quaternion Product: 



[V 


' p' 


q 2 


Pi 

q = 

<h 

p = 

Pi 


Wj 


kP4 , 


(155) 


(156) 


(157) 


q*p = 


q 2 Pi +qiP 2 ~ ( UP\ + c hP4 
<hPl+<l4P2+<llP3-<hP4 
v tflPl -(PPl + <hPl +C UP4y 


(158) 
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